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ABSTRACT 

We present a new implementation of radiation hydrodynamics (RHD) in the adaptive 
mesh refinement (AMR) code RAMSES. The multi- group radiative transfer (RT) is 
performed on the AMR grid with a first-order Godunov method using the Ml closure 
for the Eddington tensor, and is coupled to the hydrodynamics via non-equilibrium 
thermochemistry of hydrogen and helium. This moment-based approach has the large 
advantage that the computational cost is independent of the number of radiative 
sources - it can even deal with continuous regions of emission such as bound-free 
emission from gas. As it is built directly into RAMSES, the RT takes natural advantage 
of the refinement and parallelization strategies already in place. Since we use an explicit 
advection solver for the radiative transport, the time step is restricted by the speed of 
light - a severe limitation that can be alleviated using the so-called "reduced speed of 
light" approximation. We propose a rigorous framework to assess the validity of this 
approximation in various conditions encountered in cosmology and galaxy formation. 
We finally perform with our newly developed code a complete suite of RHD tests, 
comparing our results to other RHD codes. The tests demonstrate that our code 
performs very well and is ideally suited for exploring the effect of radiation on current 
scenarios of structure and galaxy formation. 
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1 INTRODUCTION 

With the surging interest in reionization and the first sources 
of light in the Universe, and also thanks to a steadily increas- 
ing computational power, cosmological simulation codes 
have begun to include ionizing radiative transfer (RT) in 
the last decade or so. This is generally seen as a second- 
order component in most astrophysical processes, but im- 
portant nonetheless, and is obviously very important in the 
context of simulating reionization. Due to the challenges 
involved, most implementations have started out with the 
post-processing of ionizing radiation on simulations includ- 
ing only dark matter, but a few have begun doing coupled 
radiation hydrodynamics (RHD), which model the interplay 
of radiation and gas. 
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It is highly desirable to follow self-consistently, with 
RHD simulations, the time-evolution and morphology of 
large-scale intergalactic medium (IGM) reionization and at 
the same time the smaller scale formation of the presumed 
sources of reionization; how galaxy formation is regulated 
by the ionizing radiation being released, how much of the 
radiation escapes from the galaxies to ionize the IGM, how 
first generation stars are formed in a metal-free environment 
and how radiative and supernovae feedback from those stars 
affect the inter-galactic medium. The galaxies and the IGM 
are indeed inter-connected via the ionizing radiation: the 
photons released from the galaxies affect the state of the 
surrounding gas via ionization and heating and may even 
prevent it from falling in or condensing into external gravi- 



tational potentials, especially small ones (e.g. Wise & Abel 



2008 Ocvirk & Aubert 2011), which can then in turn sig- 



nificantly alter the ionization history. 
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The importance of RT and RHD is of course not limited 
to the epoch of reionization. Stars keep emitting ionizing ra- 
diation after this epoch and their radiative feedback hkely 
has an effect on the post-reionization regulation of star- 
formation (e.g. Pawlik &: Schaye|2009 Hopkins, Quataert & 



MurraypQll ), the mass distribution of stellar populations 



(Krumholz, Klein & McKee 20121) and even gas outflows 
( Hopkins et al.||2Q12t . 



Radiation hydrodynamics are complex and costly in 
simulations. The inclusion of coupled radiative transfer in 
hydrodynamical codes in general is challenging mainly be- 
cause of the high dimensionality of radiative transfer (space-, 
angular-, and frequency dimensions) and the inherent differ- 
ence between the typical timescales of radiative transfer and 
non-relativistic hydrodynamics. Simulating the interaction 
between small and large scales (so relevant to the epoch of 
cosmic reionization) makes things even worse: one wants to 
simulate, in a statistically significant region of the Universe 
(i.e. of the order of 100 comoving Mpc across) the conden- 
sation of matter in galaxy groups on Mpc scales, down to 
individual galaxies on kpc scales, followed by the formations 
of stellar nurseries in those galaxies on pc scales, and ulti- 
mately the formation of stars on sub-pc scales and then the 
effect of radiation from those stars back to the large scale 
IGM. This cycle involves size differences of something like 
9 to 10 orders of magnitude - which is too much for the 
most advanced codes and computers today, actually even so 
without the inclusion of radiative transfer. 

Due to these challenges, simulations typically focus on 
only a subset of these scales; either they consider reioniza- 
tion on large scales and apply sub-resolution recipes to de- 
termine stellar luminosities and UV escape fractions, or they 
ignore the cosmological context and focus on star formation 
and escape fractions in isolated galaxies or even isolated 
stellar nurseries. 

A number of large scale 3D radiative transfer simula- 
tions of reionization have been carried out in recent years 
(e.g. iGnedin & Ostriker||1997| |Miralda-Escude, Haehnelt 



fc Rees||2000| |G nedin"2000': Ciardi , Stoehr fc White] 
Sokasian et al.|l2004i Iliev et al. 2006b| IZahn et al. 



2010 



Croft fc Altay|l2008| |Baek et al ll2010| [Aubert h Teyssier] 
Petkova fc Springel||2011b ), though they must all to 



2003 



2007 



some degree use subgrid recipes for star formation rates, stel- 
lar luminosities and UV escape fractions, none of which are 
well constrained. The ionization history in these simulations 
thus largely depends on these input parameters and resolu- 
tion - some in fact use the observational constraints of the 
ionization history to derive constraints on these free param- 



eters (e^gJSo^sian^r^ 2004; C roft fc Altay|| 2008; Back 
eFldll201Q| lAubert fc Teyssier.. 20 10| [Petkova fc Springel. 
201 lb|). Furthermore, most of these works have used a post- 



processing RT strategy instead of RHD, which neglects the 
effect the ionizing radiation has on the formation of lumi- 
nous sources. 

The primary driver behind this work is the desire to 
understand the birth of galaxies and stars during the dark 
ages, and how they link with their large scale environment. 
We have thus implemented a RHD version of the widely 
used cosmological code RAMSES ( Teyssier|2002 ), that we call 
RAMSES-RT, with the goal of running cosmological RHD sim- 
ulations, optimized for galactic scale radiation hydrodynam- 
ics. RAMSES is an adaptive mesh refinement (AMR) code, 



which greatly cuts costs by adaptively allowing the resolu- 
tion to follow the formation of structures. The RHD imple- 
mentation takes full advantage of the AMR strategy, allow- 
ing for high resolution simulations that can self consistently 
model the interplay of the reionizing Universe and the for- 
mation of the first galaxies. 

Some of the goals we will be able to tackle with this 
implementation are: 

• Study radiative feedback effects in primordial galaxies. 
These galaxies are by definition young and small, and the 
first stars are thought to be gigantic and very bright due 
to the lack of metals. The ionizing radiation from these 
first stars is likely to have a dramatic effect on the galaxy 
evolution. This is closely associated with the formation of 
molecules, needed to form the first stars, which is sensitive 
to the radiation field. Radiative feedback effects also appear 
to be relevant in lower-redshift galaxies, and likely have a 
considerable impact on the initial mass function of stellar 
populations (K rumholz, Klein fc McKee||2012j ). 

• Investigate the escape of ionizing photons from early 
galaxies, how it affects the ionization history and external 
structure formation, e.g. the formation of satellite galaxies. 

• Study the emission and absorption properties of galax- 
ies and extended structures. Observable properties of gas 
are highly dependent on its ionization state, which in turn 



depends on the local radiation field (e.g. Oppenheimer & 
Schaye '2013|. To predict it correctly, and to make correct 
interpretations of existing observations, one thus needs to 
model the ionization state consistently, for which RHD sim- 
ulations are needed. 

• Improve sub-resolution recipes: of course we have not 
implemented a miracle code, and we are still nowhere near 
simulating simultaneously the 9 to 10 orders of magnitude in 
scale needed for fully self- consistent simulations of reioniza- 
tion. Sub-resolution strategies are still needed, and part of 
the objective is to improve those via small-scale simulations 
of stellar feedback (SNe, radiation, stellar winds). 

It is useful here to make clear the distinction between 
continuum and line radiative transfer: our goal is to study 
the interplay of ionizing radiation, e.g. from stellar popula- 
tions and AGN, and the interstellar/intergalactic gas. We 
consider continuum radiation, because the spectra of stars 
(and AGN) are smooth enough that emission and absorption 
processes are not sensitive to subtle rest-frame frequency 
shifts, be they due to local gas velocities or cosmological 
expansion. 

On the other side is line transfer, i.e. the propagation 
of radiation over a narrow frequency range, usually corre- 
sponding to a central frequency that resonates with the gas 
particles. An important example is the propagation of Lya 
photons. Here, one is interested in the complex frequency 
and direction shifts that take place via scattering on the 
gas particles, and gas velocities and subtle frequency shifts 
are vital components. Line transfer is mostly done to inter- 
pret observational spectra, e.g. from Lya emitting/absorb- 
ing galaxies (e.g. Verhamme, Schaerer & Maselli"2006), and 
is usually run in post-processing under the assumption that 
the line radiation has a negligible effect on the gas dynamics 



(through this assumption is not neccessarily true; see Dijk- 
stra fc Loeb||2009t . 

There is a bit of a grey line between those two regimes 
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of continuum and line radiation - some codes are even able 
to do both (e.g.|Baek et al.|2009||Pierleoni, Maselli fc Ciardi 



2009 Yajima et al.|2012 ) . Our implementation deals strictly 
with continuum radiation though, as do most RHD imple- 
mentations, for the sake of speed and memory limitations. 
We do approximate multi-frequency, but only quite coarsely, 
such that simulated photons represent an average of photons 
over a relatively wide frequency range, and any subtle fre- 
quency shifts and velocity effects are ignored. 

1.1 Radiative transfer schemes and existing 
implementations 

Cosmological hydrodynamics codes have traditionally been 
divided into two categories: Smoothed Particle Hydrody- 
namics (SPH) and AMR. The drawbacks and advantages 
of each method have been thoroughly explored (e.g. Agertz 



et al.|2007| [Wadsley, Veeravalh fc Couchman||2508| JTasker 
et al.||2008{ and we now believe that both code types agree 
more or less on the final result if they are used carefully with 
recently developed fixes and improvements. On the radiation 
side, it is quite remarkable that we have the same dichotomy 
between ray-tracing codes and moment-based codes. Com- 
parative evaluations of both methods have been performed 



in several papers (Iliev et al. 2006a Alt ay, Croft & Pelu- 



pessy|20 08; Aubert & Teyssier 2008; Pawlik & Schaye 2008; 
Iliev et al..2009,.Maseni, Ciardi fc Kanekar,2009, Petkova fc 



Springel|2009||Cantalupo fc Porciani|2011|(Pawlik fc Schaye 



2011||Petkova fc Springel||2011a| [Wise fc Abel||2011t , and 
here again, each method has its own specific advantage over 
the other one. Comparing both methods in the coupled case 
(RHD) within the more challenging context of galaxy for- 
mation, such as in the recent Aquila comparison project, 
(Scannapieco et aL]|2012), remains to be done. 



1.1.1 Ray-based schemes 

Here the approximation is made that the radiation field is 
dominated by a limited number of sources. This allows one to 
approximate the local intensity of radiation, /^y, as a function 
of the optical depth r along rays from each source. 

The simplest solution is to cast rays, or long charac- 
teristics from each source to each cell (or volume element) 
and sum up the optical depth at each endpoint. With the op- 
tical depths in hand, lu is known everywhere and the rates 
of photoionization, heating and cooling can be calculated. 
While this strategy has the advantage of being simple and 
easy to parallelize (each source calculation is independent 
from the other) , there is a lot of redundancy, since any cell 
which is close to a radiative source is traversed by many rays 
cast to further-lying cells, and is thus queried many times 
for its contribution to the optical depth. The parallelization 
is also not really so advantageous in the case of multipro- 
cessor codes, since rays that travel over large lengths likely 
need to access cell states over many CPU nodes, calling for a 
lot of inter-node communication. Furthermore, the method 
is expensive: the computational cost scales linearly with the 
number of radiative sources, and each RT timestep has order 
C'(A^sources A^ceiis) opcratlous, whcrc A^sources is the number 
of radiative sources and A^ceiis is the number of volume el- 
ements. Implementation examples in clude [Abel, Norman &:] 
Madau| ( [l999| ), |Cen| ( |2002| and [Susa| ( [2006| ^ 



Short characteristics schemes overcome the redun- 
dancy problem by not casting separate rays for each desti- 
nation cell. Instead, the calculation of optical depths in cells 
is propagated outwards from the source, and is in each cell 
based on the entering optical depths from the inner-lying 
cells. Calculation of the optical depth in a cell thus requires 
some sort of interpolation from the inner ones. There is no 
redundancy, as only a single ray segment is cast through 
each cell in one time-step. However, there is still a large 
number of operations and the problem has been made in- 
herently serial, since the optical depths must be calculated 
in a sequence which follows the radiation ripple away from 
the source. Some examples are|Nakamoto, Umemura fc Susa| 



Alvarez, Bromm & Shapiro (2006) 



eJN^ 
jWh 



([200T|),periema et al.|(p006 ) [Whalen fc Norman, ( ,2006J and 



Adaptive ray tracing (e.g. Abel & Wandelt 2002 



Razoumov fc Cardall||2005| [Wise fc Abel||2011t is avariant 
on long characteristics, where rays of photons are integrated 
outwards from the source, updating the ray at every step of 
the way via absorption. To minimize redundancy, only a 
handful of rays are cast from the source, but they are split 
into sub-rays to ensure that all cells are covered by them, 
and they can be merged again if need be. 

Cones are a variant on short characteristics, used in 
conjunction with SPH ( jPawlik fc Sch aye 20 08| |2011[) and 
the moving- mesh AREPO code (Petk ova fc Springel|2011a" ). 
The angular dimension of the RT equation is discretized into 
tesselating cones that can collect radiation from multiple 
sources and thus ease the computational load and even allow 
for the inclusion of continuous sources, e.g. gas coUisional 
recombination. 

A hybrid method proposed bylRijkhorst et al. (2006) 



combines the long and short characteristics on patch-based 
grids (like AMR), to get rid of most of the redundancy while 
keeping the parallel nature. Long characteristics are used 
inside patches, while short characteristics are used for the 
inter-patches calculations. 

Monte- Carlo schemes do without splitting or merging 
of rays, but instead reduce the computational cost by sam- 
pling the radiation field, typically both in the angular and 
frequency dimensions, into photon packets that are emit- 
ted and traced away from the source. The cost can thus 
be adjusted with the number of packets emitted, but gen- 
erally this number must be high in order to minimize the 
noise inherent to such a statistical method. Examples in- 



clude Ciardi et al. (2001), Maselh, Ferrara k Ciardi (2003), 



Altay, Croft fc Pelupessy] (|2008| ), |Baek et al.| ( |2009t , and 
Cantalupo fc Porciani| ( |2011| ). An advantage of the Monte- 
Carlo approach of tracking individual photon packets is that 
it naturally allows for keeping track of the scattering of pho- 
tons. For line radiation transfer, where doppler/redshift ef- 
fects in resonant photon scattering are important, Monte- 
Carlo schemes are the only feasible way to go - though 
in these cases, post-processing RT is usually sufficient (e.g. 



Cantalupo et al.|2005||Verhamme, Schaerer &: Maseni|2006| 
Laursen fc Sommer-Larsen|2007||Pierleoni, Maselli fc Ciardi| 



2009| ). 



Ray-based schemes in general assume infinite light 
speed, i.e. rays are cast from source to destination instanta- 
neously. Many authors note that this only affects the initial 
speed of ionization fronts (I-fronts) around points sources 
(being faster than the light speed), but it may also result 



© 0000 RAS, MNRAS 000, 000-000 



Rosdahl et al. 



in an over-estimated I-front speed in underdense regions 



(see ^6.5), and may thus give incorrect results in reion- 



ization experiments where voids are re-ionized too quickly. 
Some ray schemes (e.g. Pawlik fc Schaye 2008 



Springel 2011a Wise k Abel 2011) allow for finite light 



Petkova & 



speed, but this adds to the complexity, memory require- 
ment and computational load. It is also unavoidable with 
ray-based schemes that the computational load increases 
with the number of radiative sources. This problem largely 
disappears with moment methods, though of course other 
limitations appear instead. 

1.1.2 Moment-based radiative transfer 

An alternative to ray-tracing schemes is to reduce the angu- 
lar dimensions by taking angular moments of the RT equa- 
tion (Eq.|2]). Intuitively this can be thought of as switching 
from a beam description to that of a field or a fiuid, where 
the individual beams are replaced with a "bulk" direction 
that represents an average of all the photons crossing a given 
volume element in space. This infers useful simplifications: 
two angular dimensions are eliminated from the problem, 
and the equations take a form of conservation laws which 
are akin to the Euler equations of hydrodynamics, and are 
thus rather easily coupled to these equations, and can even 
be solved with numerical methods designed for hydrodynam- 
ics. 

The main advantage is also the main drawback: the di- 
rectionality is largely lost in the moment approximation and 
the radiation becomes somewhat diffusive, which is gener- 
ally a good description in the optically thick limit, where 
the radiation scatters a lot, but not in the optically thin 
regime where the radiation is free- streaming. Radiation has 
a tendency to creep around corners with moment methods. 
Shadows are usually only coarsely approximated, if at all, 
though we will see e.g. in section §6.4| that sharp shadows 
can be maintained with idealized setups and a specific solver. 

The large value of the speed of light is also an issue. Mo- 
ment methods based on an explicit time marching scheme 
have to follow a Courant stability condition that basically 
limits the radiation from crossing more than one volume ele- 
ment in one time-step. This requires to perform many time- 
steps to simulate a light crossing time in the free- streaming 
limit, or, as we will see later, to reduce artificially the light 
speed. Implicit solvers can somewhat alleviate this limita- 
tion, at the price of inverting large sparse matrices which are 
usually ill-conditioned and require expensive, poorly paral- 
lelized, relaxation methods. 

In moment based implementations, the frequency di- 
mension is also reduced, via integration over frequency bins: 
in the grey (single group) approximation the integral is per- 
formed over the whole relevant frequency range, typically 
from the hydrogen ionization frequency and upwards. In the 
multigroup approximation, the frequency range is split into 
a handful of bins, or photon groups, (rarely more than a few 
tens due to memory and computational limitations) and the 
equations of radiative transfer can be solved separately for 
each group. 

In the simplest form of moment-based RT implementa- 
tions, so-called fiux limited diffusion (FLD), only the zero- 
th order moment of the radiative transfer equation is used, 
and a closure is provided in the form of a local diffusion 



relation, which lets the radiation fiow in the direction of de- 
creasing gas internal energy (i.e. in the direction opposite 
of the energy gradient). This is realistic only if the medium 
is optically thick, and shadows cannot be modelled. The 



FLD method has been used by e.g. Krumholz et al. (2007), 
Reynolds et al. (2009') and'Commergon et al. (2011 ), mainly 



for the purpose of studying the momentum feedback of in- 
frared radiation onto dusty and optically thick gas, rather 
than photoionization of hydrogen and helium. 



Gnedin fc Abel (2001 ) and Petkova k Springel (2009) 



used the optically thin variable Eddington tensor formalism 
(OTVET), in which the direction of the radiative field is 
composed on-the-fiy in every point in space from all the ra- 
diative sources in the simulation, assuming that the medium 
between source and destination is transparent (hence opti- 
cally thin). This calculation is pretty fast, given the num- 
ber of relevant radiative sources is not ov erburdening, and 
one can neglect these in-between gas cells. Finlator, Ozel & 



Dave (|2009|) take this further and include in the calculation 
the optical thickness between source and destination with 
a long characteristics method, which makes for an accurate 
but slow implementation. 



Gonzalez, Audit k Huynh (2007), AT08 and Vaytet, 



Audit & Dubroca (2010) - and now us - use a different 



closure formalism, the so-called Ml closure, which can es- 
tablish and retain bulk directionality of photon fiows, and 
can to some degree model shadows behind opaque obstacles. 
The Ml closure is very advantageous in the sense that it is 
purely local, i.e. it requires no information which lies outside 
the cell, which is not the case for the OTVET approxima- 
tion. 



As shown by Dubroca & Feugeas (1999), the Ml clo- 
sure has the further advantage that it makes the system of 
RT equations take locally the form of a hyperbolic system 
of conservations laws, where the characteristic wave speeds 
can be calculated explicitly and are usually close, but al- 
ways smaller than the speed of light c. Hyperbolic systems 
of conservation laws are mathematically well understood and 
thoroughly investigated, and a plethora of numerical meth- 
ods exist to deal with them (e.g. |Toro||1999|. In fact, the 
Euler equations are also a hyperbolic system of conserva- 
tion laws, which implies we have the RT equations in a form 
which is well suited to lie alongside existing hydrodynamical 
solvers, e.g. in RAMSES. 

1.2 From ATON to RAMSES-RT 



The ATON code (AT08, |Aubert fc TeyssierpOlOl uses graphi- 
cal processing units, or GPUs, to post-process the transfer of 
monochromatic photons and their interaction with hydrogen 
gas. GPUs are very fast, and therefore offer the possibility to 
use the correct (very large) value for the speed of light and 
perform hundreds to thousands of radiation sub-cycles at a 
reasonable cost, but only if the data is optimally structured 
in memory, such that volume elements that are close in space 
are also close in memory. It is ideal for post-processing RT 
on simulation outputs that are projected onto a Cartesian 
grid, but hard to couple directly with an AMR grid in or- 
der to play an active part in any complex galaxy formation 
simulation. Even so, we have in the newest version of the 
ATON code included the possibility to perform fully coupled 
RHD simulations using a Cartesian grid only (this usually 
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corresponds to our coarser grid level in the AMR hierarchy) , 
where RT is performed using the AXON module on GPUs. 

In our RAMSES-RT implementation we use the same RT 
method as ATON does - the moment method with the Ml 
Eddington tensor closure. The biggest difference is that 
RAMSES-RT is built directly into the RAMSES cosmological hy- 
drodynamics code, allowing us to perform RHD simulations 
directly on the AMR grid, without any transfer of data be- 
tween different grid structures. Furthermore, we have ex- 
panded the implementation to include multigroup photons 
to approximate multifrequency, and we have added the in- 
teractions between photons and helium. We explicitly store 
and advect the ionization states of hydrogen and helium, 
and we have built into RAMSES-RT a new non-equilibrium 
thermochemistry model that evolves these states along with 
the temperature and the radiation field through chemical 
processes, photon absorption and emission. Finally, for re- 
alistic radiative feedback from stellar populations, we have 
enabled RAMSES-RT to read external spectral energy distri- 
bution (SED) models and derive from them luminosities and 
UV "colors" of simulated stellar sources. 

We have already listed a number of RT implementa- 
tions, two of which even function already in the RAMSES 



code (AT08, Commergon et al. 2011), and one might ask 
whether another one is really needed. To first answer for the 
ATON implementation, it is optimized for a different regime 
than RAMSES-RT. As discussed, ATON prefers to work with 
structured grids, but it cannot deal well with adaptive re- 
finement. This, plus the speed of ATON, makes it very good 
for studying large scale cosmological reionization, but not 
good for AMR simulations of individual halos/galaxies, e.g. 
cosmological zoom simulations, where the subject of interest 
is the effect of radiative feedback on the formation of struc- 
tures and galaxy evolution, and escape fractions of ionizing 



radiation. The Commergon et al. (2011) implementation is 



on the opposite side of the spectrum. Being based on the 
FLD method, it is optimized for RHD simulations of op- 
tically thick protostellar gas. It is a monogroup code that 
doesn't track the ionization state of the gas. Furthermore, 
it uses a rather costly implicit solver, which makes it hard 
to adapt to multiple adaptive time stepping usually used in 
galaxy formation problems. 

Out of the other codes we have listed, only three of 
those seem to have been used for published 3D cosmologi- 
cal RHD simulations with ionizing radiation. As far as we 
can see these are Gnedin & Abel (2001) (in ART) , [Petkova 



fc Springell ( [2Q09t (in GADGET), and [Wise &; Abel| ( p01lf (in 
ENZO). A few others that have been used for published as- 
trophysical (ionizing) RHD simulations but without a co- 
evolving cosmology are Mellema et al. ( 2006 ) , Susa ( 2006 ) , 



Whalen & Norman (2006), and Baek et al. (2009). The rest 



apparently only do post-processing RT, aren't parallel or are 
otherwise not efficient enough. Many of these codes are also 
optimized for cosmological reionization rather than galaxy- 
scale feedback. 

Thus there aren't so many cosmological RHD imple- 
mentations out there, and there should be room for more. 
The main advantage of our implementation is that our 
method allows for an unlimited number of radiative sources 
and can even easily handle continuous sources, and is thus 
ideal for modelling e.g. the effects of radiative feedback in 
highly resolved simulations of galaxy formation, UV escape 



fractions, and the effects of self-shielding on the emission 
properties of gas and structure formation, e.g. in the con- 
text of galaxy formation in weak gravitational potentials. 

The structure of the paper is as follows: in ^we present 
the moment based RT method we use. In ^we explain how 
we inject and transport radiation on a grid of gas cells, and 
how we calculate the thermochemistry in each cell, that in- 
corporates the absorption and emission of radiation. In El 
we present two tricks we use to speed up the RHD code, 
namely to reduce the speed of light, and to "smooth" out 
the effect of operator splitting. In ^we describe how the ra- 
diative transfer calculation is placed in the numerical scheme 
of RAMSES, and demonstrate that the radiation is accurately 
transported across an AMR grid. In ^ we present our test 
suite, demonstrating that our code performs very well in 
coupled radiation hydrodynamics problems and finally, ^ 
summarizes this work and points towards features that may 
be added in the future. Details of the thermochemistry and 
additional code tests are described in the appendix. 



2 MOMENT-BASED RADIATIVE TRANSFER 
WITH THE Ml CLOSURE 

Let /iy(x, n,t) denote the radiation specific intensity at lo- 
cation X and time t, such that 



/jy dv dQ dA dt 



(1) 



is the energy of photons with frequency over the range du 
around u propagating through the area dA in a solid angle 
dQ around the direction n. 

The equation of radiative transfer (e.g. Mihalas fc Mi- 
halas|1984 ) describes the local change in lu as a function of 
propagation, absorption and emission. 



c dt 



+ n • VIu = -K,ulu + r]u, 



(2) 



where c is the speed of light, K,u{^,n,t) is an absorption 
coefficient and ?7i/(x, n, t) a source function. 

By taking the zeroth and first angular moments of ([2|, 
we can derive the moment-based RT equations that describe 
the time-evolution of photon number density N^ and flux Fj^ 
(see e.g. AT08): 



dt 



Hl,Hel,Hell 



+ V-F.. = - Yl rija^^jcN^ + Nt + K^"" (3) 



dt 



Hl,HeI,HeII 



+ c^V • P^ = - ^ rija^jcY^, 



(4) 



where Pj/ is the radiative pressure tensor that remains to 
be determined to close the set of equations. Here we have 
split the absorption coefficient into constituent terms, nja^j, 
where rij is number density of the photo-absorbing species j 
(=Hi, Hei, Hell), and ai,j is the ionization cross section be- 
tween z/- frequency photons and species j. Furthermore we 
have split the source function into (e.g. stellar, quasar) in- 
jection sources, N^^, and recombination radiation from gas, 
N^^^. Here we only consider the photo-absorption of hy- 
drogen and helium, which is obviously most relevant in 
the regime of UV photons. However, other absorbers can 
straightforwardly be added to the system. 
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Eqs. ([3|-(|4| are continuous in zy, and they must be dis- 
cretized to be usable in a numerical code. AT08 collected all 
relevant frequencies into one bin, so the equations could be 
solved for one group of photons whose attributes represent 
averages over the frequency range. For a rough approxima- 
tion of multifrequency, we split the relevant frequency range 
into a number of photon groups, defined by 






Nr, dv, 






diy, 



(5) 



where (i^io,^ii) is the frequency interval for group i. In the 
limit of one photon group, the frequency range is {fiQ^fn) — 
(z/Hi, oo); with M > 1 groups, the frequency intervals should 
typically be mutually exclusive and set up to cover the whole 
H-ionizing range: 

[z^oo,z^oi : z^io,z^ii : ••• : t^mo.t^mi] — [z^hi,oo[. 

Integrating the RT equations ([3| and (El over each fre- 
quency bin corresponding to the group definitions yields M 
sets of four equations: 



dt 



Hl,Hel,Hell 



+ V.F. = - Yl rijCafjN, + iV^ + iV^^ (6) 



Hl,HeI,HeII 
oFi 2t-7 -IT) V^ N-c^ 

— -+cV-Pi = - 2^ UjCa.jFi 



(7) 



where afj represent average cross sections between each 



group i and species j, defined b}jj 

p^ N, dv ' 



(8) 



We simplify things however by defining the group cross sec- 
tions as global quantities, assuming a frequency distribution 
of energy J{y) for the radiative sources (e.g. a black-body 
or some sophisticated model) . The cross sections are thus in 
practice evaluated by 






i::: Ji^)/^^ ^^ 



(9) 



where hv is photon energy (with h the Planck constant). 
Likewise, average photon energies within each group are 
evaluated by 



ei 



/;^;^ j{v) dv 

s:{j{v)/hvdv' 



(10) 



and furthermore, for the calculation of photoionization heat- 
inqj energy weighted cross sections are stored for each group 
- absorbing species couple: 



n^ J{iy) diy ■ 



(11) 



In RAMSES-RT, afj, afj and Ci can be either set by hand or 
evaluated on-the-fly from spectral energy distribution tables 
as luminosity weighted averages from in-simulation stellar 



^ here we assume the spectral shape of Fjy to be identical, within 
each grou p, to that of N^. 
2 see Eq. (a16| 



populations, using the expressions from Verner et al. (1996) 
for cri.,Hi, cru,uei and au^ueii- 

For each photon group, the corresponding set of equa- 
tions ([6|-([7|) must be closed with an expression for the pres- 
sure tensor P. This tensor is usually described as the prod- 
uct of the photon number density and the so-called Edding- 
ton tensor D (see Eq. 12), for which some meaningful and 
physical expressio n is desired . Some formalisms have been 
suggested for D^- 



Dave 



Gnedin k Abel 



(2001), 



(2009), and Petkova fc Springel ( 2009 ) have used the so 



Finlator, Ozel & 



called optically thin Eddington tensor formalism (OTVET), 
in which P is composed on-the-fly from all the radiation 
sources, the main drawback being the computational cost 
associated with collecting the positions of every radiative 
source relative to every volume element. Instead, like AT08 
(and Gonzalez, Audit fc Huynh||2007 [ before them) , we use 
the Ml closure relation (|Levermore| 1984), which has the 
great advantages that it is purely local, i.e. evaluating it in 
a piece of space only requires local quantities, and that it 
can retain a directionality along the flow of the radiative 
field. In our frequency-discretized form, the pressure tensor 
is given in each volume element for each photon group by 



P^ 



■■Tl^N^ 



where the Eddington tensor is 



D, 



Xi 



1 + 



3Xi 



and 



ni = 



m 



Xi = 



3 + 4/^ 



5 + 2^4^3^ 



f^ = 



CN^ 



(12) 



(13) 



(14) 



are the unit vector pointing in the flux direction, the Edding- 
ton factor and the reduced flux, respectively. The reduced 
flux describes the directionality of the group i radiation in 
each point, and must always have ^ /^ ^ 1. A low value 
means the radiation is predominantly isotropic, and a high 
value means it is predominantly flowing in one direction. For 
the arguments leading to these expressions and a general dis- 



cussion we point the reader towards Levermore (1984) and 
Gonzalez, Audit k Huynh| ( |2007| and AT08. 



3 THE RADIATIVE TRANSFER 
IMPLEMENTATION 

We will now describe how pure radiative transfer is solved 
on a grid - without yet taking into consideration the hydro- 
dynamical coupling. The details here are not very specific 
to RAMSES-RT and are much hke those of AT08. 

In addition to the usual hydrodynamical variables 
stored in every grid cell in RAMSES (gas density p, momen- 
tum density pu, energy density E, metallicity Z), RAMSES-RT 
has the following variables: First, we have the 4 x M vari- 
ables describing photon densities Ni and fluxes F^ for the M 
photon groups. Second, in order to consistently treat the in- 
teractions of photons and gas, we track the non-equilibrium 
evolution of hydrogen and helium ionization in every cell, 
stored in the form of passive scalars which are advected with 
the gas, namely 

a^HlI = , a^Hell = , a^Helll = • (15) 
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For each photon group, we solve the set of equations 
(I6|-(l7| with an operator sphtting strategy, which involves 
decomposing the equations into three steps that are exe- 
cuted in sequence over the same time-step At, which has 
some pre-determined length. The steps are: 

(i) Photon injection step, where radiation from stellar and 
other radiative sources (other than gas recombinations) is 
injected into the grid. This corresponds to the Ni term in 

§. 

(ii) Photon transport step, where photons are propagated 
in space. This corresponds to solving ([6|-([7| with the RHS 
being equal to zero. 

(iii) Thermochemistry step, where the rest of the RHS 
of (|6}-([7| is solved. This is where the photons and the gas 
couple, so here we evolve not only the photon densities and 
fluxes, but also the ionization state and temperature of the 
gas. 



3.1 The injection step 

The equations to solve in this step are very simple. 



ON, 
dt 



--n: 



(16) 



where Ni is a rate of photon injection into photon group 
i, in the given cell. Normally, the injected photons come 
from stellar sources, but they could also include other point 
sources such as AGN, and also pre-defined point sources or 
even continuous "volume" sourceqf] 

Given the time t and time-step length At, the discrete 
update in each cell done for each photon group is the fol- 
lowing sum over all stellar particles situated in the cell: 



Nr 



Nr 



V 



(17) 



^ m. [^.(r.-+^Z.)-^.(r.^Z.)], 



where n denotes the time index {n — t and n + 1 = t + At) , 
/esc is an escape fraction, V is the cell volume, m^, r^ and Z^ 
are mass, age and metallicity of the stellar particles, respec- 
tively, and Hi is some model for the accumulated number of 
group i photons emitted per solar mass over the lifetime (so 
far) of a stellar particle. The escape fraction, /esc is just a 
parameter that can be used to express the suppression (or 
even boosting) of radiation from processes that are unre- 
solved inside the gas cell. 

RAMSES-RT can read stellar energy distribution (SED) 
model tables to do on-the-fiy evaluation of the stellar parti- 
cle luminosities, 11^. Photon cross sections and energies can 
also be determined on-the-fiy from the same tables. Details 
are given in Appendix [B] 



^ In |Rosdahl fc BlaizQt|( |2012| , we emitted UV background radi- 
ation from cosmological void regions, under the assumption that 
they are transparent to the radiation. 



3.2 The transport step 

The equations describing free-fiowing photons are 



dN 
~dt 

OF 

'dt 



+ V-F = 0, 



+ cV' 



0, 



(18) 
(19) 



i.e. (|6|-([7| with the RHS = 0. Note that we have removed 
the photon group subscript, since this set of equations is 
solved independently for each group over the time- step. 
We can write the above equations in vector form 



f+V^(«) 



0, 



(20) 



where U = [A^, F] and J^{U) = [F,c^P]. To solve ([20| over 
time-step At, we use an explicit conservative formulation, 
expressed here in ID for simplicity. 



At 



+ 



'M + 1/2 



■ 'M-1/2 



Ax 



= 0, 



(21) 



where n again denotes time index and / denotes cell index 
along the x-axis. J^i^i/2 and J-'i-i/2 = J~'(i-i)-\-i/2 are inter- 
cell fiuxes evaluated at the cell interfaces. Simple algebra 
gives us the updated cell state, 

^"+i=wr + |^(^_i/2-^+i/2), (22) 

and all we have to do is determine expressions for the inter- 
cell fiuxes. 

Many intercell fiux functions are available for differen- 
tial equations of the form (|20| which give stable results in 



bro|l{ 



the form of (22) (see e.g. Toro|1999 ), as long as the Courant 
time-step condition is respected (see {4.1). Following AT08 
and 'Gonzalez, Audit & Huynh ( 2007| we implement two fiux 
functions which can be used in RAMSES-RT. 

One is the Harten-Lax-van Leer (HLL) fiux function 
(Harten, Lax &: Leer||1983), 



(J^H 



A+ JT - A- j-r+1 + A+A- (^r+i - K) 



Jl-\-l/2 



(23) 



where 



A+ -A- 

A+ /rv A max \niax\ 

^ = max(0,Ai , A^+i ), 

A— ' { r\ A min \niin\ 

A =mm(0,Ai ,\i^x) 

are maximum and minimum eigenvalues of the Jacobian 
dT jdhi. These eigenvalues mathematically correspond to 
wave speeds, which in the case of 3D radiative transfer de- 



pend only on the magnitude of the reduced fiux / ( 14 ) and 



the angle of incidence of the fiux vector to the cell inter- 
face. This dependence has been calculated and tabulated by 
Gonzalez, Audit h Huynh (2007), and we use their table to 
extract the eigenvalues. 

The other fiux function we have implemented is the sim- 
pler Global Lax Friedrich (GLF) function. 



(J^G 






)l + l/2 



{ur+i-un 



(24) 



which corresponds to setting the HLL eigenvalues to the 
speed of light, i.e. A~ = — c and A^ = c, and has the effect 
of making the radiative transport more diffusive. Beams and 
shadows are therefore better modelled with the HLL fiux 
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Figure 1. Comparison of radiative transport in 2D RAMSES-RT runs (no gas interaction) with the HLL (top) and GLF (bottom) flux 
functions, using isotropic point sources and beams. The box width in all runs is 1 cm and the resolution is 256^ cells. For each isotropic 
point source, 10-*^° photons s"-*^ are injected continuously, and for each beam a constant photon density oi N = 1 cm~^, with a unity 
reduced flux, is imposed on a region of one cell width and eight cell heights at the beam origin. The snapshots are taken after a few 
light crossing times. Far left frames show single isotropic point sources. Middle left frames show attempts at creating horizontal and 
diagonal beams (with F = (cA/", 0) and F = {cN, cA^)/\/2, respectively). Middle right frames show two isotropic point sources and how 
the photons behave between them. Far right frames show two beams of opposing directions and how a spurious weak perpendicular 
radiation source forms where they meet. 



function than with the GLF one, whereas the inherent di- 
rectionality in the HLL function results in radiation around 
isotropic sources (e.g. stars) which is noticeably asymmetric, 
due to the preference of the axis directions. 

Fig. [l] illustrates the difference between the two flux 
functions in some idealized 2D RAMSES-RT tests, where we 
shoot off beams and turn on isotropic sources. It can be seen 
that the HLL flux function fails to give isotropic radiation 
(far left) and that the GLF function gives more diffusive 
beams (second from left). Note also how the diffusivity of 
beams with the HLL flux function is direction-dependent. 
A horizontal or vertical beam is perfectly retained while a 
diagonal one "leaks" to the sides almost as much as with 
the GLF function, which has the advantage of being fairly 
consistent on whether the beam is along-axis or diagonal. 
The right frames of the figure give an idea of how the ra- 
diative transport behaves in the case of multiple sources, 
i.e. with opposing beams and neighboring isotropic sources. 
The two opposing beams example is a typical configura- 
tion where the Ml closure relation obviously fails, creating 
a spurious source of radiation, perpendicular to the beam 
direction. However, this configuration is not really relevant 
in any astrophysical contexts, so this doesn't present a se- 
rious problem in our implementation. We generally prefer 
to use the GLF flux function, since we mostly deal with 
isotropic sources in our cosmological/galactic simulations, 
but the choice of function really depends on the problem. 
There is no noticeable difference in the computational load, 
so if shadows are important, one should go for HLL. AT08 
have compared the two flux functions in some of the bench- 



very similar results. We do likewise for the test we describe 
in §6.71 and come to the same conclusion. 



3.3 The thermochemical step 

Here we solve for the interaction between photons and gas. 
This is done by solving ([6| and ([7| with zero divergence and 
stellar injection terms. 

Photon absorption and emission have the effect of heat- 
ing and cooling the gas, so in order to self-consistently im- 
plement these interactions, we evolve along with them the 
thermal energy density e of the gas and the abundances of 
the species that interact with the photons, here Hi, Hei and 
Hell via photoionizations and Hii, Hell (again) and Hem via 
recombinations. We follow these abundances in the form of 
the three ionization fractions xhii, a^Heii and XHeiii, that we 
presented in Eqs. (15). The set of non-equilibrium thermo- 



chemistry equations solved in RAMSES-RT consists of: 



dt 



Hl,Hel,Hell 



^ TijCcr^jNi 



(25) 



Hii, Hell, Hem 



E7 T&c r A B 1 



mark RT tests of Iliev et al. ( 2006a ) and found that they give 



dt 
dt 



Hl,Hel,Hell 



y^ UjCcr^jYi, 



n + c 



(26) 
(27) 
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riH- 



dxY 



dt 



( "^ N ^ 



■ TT'HII ^HII ^e, 



r^He ^^ = ^Hel I /^Hel ^e + ^ CT^HelCA^i I 



(28) 



(29) 



+ riHelll «HeIII ^e 



riHe 



^Xt 



at 



riHell I /^Helirie + ^Hell^e + ^ CTmellCA^i I 

/ "^ N ^ 

ell ( /^Hell ^e + ^ CffllellCNi j (^ 



riH 



^Helll Q^Helll ^e 



In the photon density and flux equations, (25) and 
(26|), we have replaced the photon emission rate Nl^^ with 



the full expression for recombinative emissions from gas. 
Here, aj (T) and (y^(T) represent case A and B recom- 
bination rates for electrons combining with species j {— 
Hii, Hell, Hem). The ^Jf" factor is a boolean (1/0) that states 
which photon group j-species recombinations emit into, and 
Tie is electron number density (a direct function of the H 
and He ionization states, neglecting the contribution from 
metals) . 

The temperature-evolution, (27) is greatly sim- 
plified here (see Appendix [A] for details). Basi- 
cally it consists of two terms: the photoheating rate 
'H(A/'i,XHii,XHeii,a:Heiii,''T'H) and the radiative cooling rate 

£(T, Ni, XHII, XHell, XHelll, ^h)- 



The xhii evolution ( 28 ) consists of, respectively on the 
RHS, Hi collisional ionizations. Hi photo- ionizations, and 
Hll recombinations. Here, /3(T) is a rate of collisional ion- 



izations. The XHeii evolution (29) consists of, from left to 
right, Hel collisional ionizations. Hem recombinations, Hei 
photo-ionizations, and Hell collisional ionizations, recombi- 
nations, and photoionizations. Likewise, the XHeiii evolution 
( 30 ) consists of Hell collisional ionizations and photoioniza- 
tions, and Hem recombinations. The expressions we use for 
rates of recombinations and collisional ionizations are given 
in Appendix [E] 

The computational approach we use to solving Equa- 
tions (25)- (30) takes inspiration from Anninos et al. (1997). 



The basic premise is to solve the equations over a sub-step in 
a specific order (the order we have given) , explicitly for those 
variables that remain to be solved (including the current 
one), but implicitly for those that have already be solved 



over the sub-step. Eqs. (25) and (26) are thus solved purely 



explicitly, using the backwards-in-time (BW) values for all 
variables on the RHS. Eq. (27) is partly implicit in the sense 
that it uses forward-in-time values for N and F, but BW 
values for the other variables. And so on, ending with Eq. 
( 30 ) , which is then implicit in every variable except the one 
solved for (xHeiii)- We give details of the discretization of 
these equations in Appendix [A] 



3.3.1 The 10% thermochemistry rule 

For accuracy, each thermochemistry step is restricted by a 
local cooling time which prohibits any of the thermochem- 
ical quantities to change by a substantial fraction in one 



time-step. We therefore sub-cycle the thermochemistry step 
to fill in the global RT time-step (see next section), using 
what can be called the 10% rule: In each cell, the thermo- 
chemistry step is initially executed with the full RT time- 
step length, and then the fractional update is considered. If 
any of the evolved quantities (A/'i, F^, e, ionization fractions) 
have changed by more than 10%, we backtrack and do the 
same calculation with half the time-step length. Conversely, 
if the greatest fractional change in a sub-step is < 5%, the 
timestep length is doubled for the next sub-step (without 
the backtracking). 



3.3.2 The on-the-spot approxmation 

The photon-emitting recombinative term, the second RHS 



sum in ( 25 ), is optionally included. Excluding it is usually re- 



ferred to as the on-the-spot approximation (OTSA), mean- 
ing that any recombination-emitted photons are absorbed 
"on the spot" by a near- lying atom (in the same cell), and 
hence these photon emissions cancel out by local photon ab- 
sorptions. If the OTSA is assumed, the gas is thus not pho- 
toemitting, and the case A recombination rates are replaced 
with case B recombination rates in (25)- (30), i.e. photon- 



emitting recombinations straight down to the ground level 
are not counted. The OTSA is in general a valid approxi- 
mation in the optically thick regime but not so when the 
photon mean free path becomes longer than the cell width. 
It is a great advantage of our RT implementation (and 
of moment-based RT in general) that it is not restricted 
to a limited number of point sources. The computational 
load does not scale at all with the number of sources, and 
photon emission from gas (non-OTSA) comes at no added 
cost, whereas it may become prohibitively expensive in ray- 
tracing implementations. 



4 TIMESTEPPING ISSUES 

RT is computationally expensive, and we use two basic tricks 
to speed up the calculation. One is to reduce the speed of 
light, the other is to modify slightly the traditional opera- 
tor splitting approach, by increasing the coupling between 
photon injection and advection on one hand and thermo- 
chemistry and photo- heating on the other hand. 



4.1 The RT timestep and the reduced speed of 
light 

In each iteration before the three RT steps of photon in- 
jection, advection and thermochemistry are executed, the 
length of the time-step, AtRT, must be determined. 

We use an explicit solver for the radiative transport 
(21), so the advection timestep, and thus the global RT 



timestep, is constrained by the Courant condition (here in 
3D), 

Ax 



AtRT < 



3c 



(31) 



where Ax is the cell width. This time-step constraint is se- 
vere: it results in an integration step which is typically 300 
times shorter than in non-relativistic hydrodynamical simu- 
lations, where the speed of light is replaced by a maximum 
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gas velocity (^ 1000 km/ s) in Eq. 31 In a coupled (RHD) 
simulation, this would imply a comparable increase in CPU 
time, either because of a global timestep reduction (as we 
chose to implement, see Sec. [5|, or because of many radia- 
tive sub-steps (as is implemented e.g. in ATOrfj). In the case 
of radiative transfer with the moment equations, there are 
two well-known solutions to this problem. 

The first solution is to use an implicit method instead 
of an explicit one to solve the transport equation, which 
means using forward-in-time intercell fluxes in (21), i.e re- 
placing F"^ = T^ with T"^^^ = J'*+^^ This seemingly sim- 
ple change ensures that the computation is always stable^ 
no matter how big the time- step, and we can get rid of 
the Courant condition. However: (i) It doesn't mean that 
the computation is accurate, and in fact we still need some 
time-stepping condition to retain the accuracy, e.g. to re- 
strain any quantity to be changed by more than say 10% 
in a single time-step. Furthermore, such a condition usually 
must be checked by trial-and-error, i.e. one guesses a time- 
step and performs a global transport step (over the grid) 
and then checks whether the accuracy constraint was broken 
anywhere. Such trial-and-error time-stepping can be very ex- 
pensive since it is a global process, (ii) Replacing T^ with 
■pt+At -g g^ctually not simple at all. Eqs. 21 become a sys- 
tem of coupled algebraic equations that must be solved via 
matrix manipulation in an iterative process, which is com- 
plicated, computationally expensive, and of limited scope 
(i.e. can't be easily applied to any problem). Due to these 
two reasons we have opted out of the implicit approach. It 
is absolutely a valid approach however, and used by many 
(e.g. Petkova & Springel 2009; Comme rgon et al^|2011 ). 

The second solution, which we have chosen instead, is 
to keep our solver explicit, and relax the Courant condi- 
tion by changing the speed of light to a reduced light speed 
Or <C c, the payoff being that the time- step (31) becomes 



longer. This is generally referred to as the reduced speed of 



light approximation (RSLA), and was introduced by ,Gnedin 
k Abel|([200T]). The idea of the RSLA is that in many appli- 



cations of interest, the propagation of light is in fact limited 
by the much slower speed of ionizing fronts. In such situ- 
ations, reducing the speed of light, while keeping it higher 
than the fastest I- front, will yield the correct solution at a 
much reduced CPU cost. In the following section, we provide 
a framework to help judge how accurate the RSLA may be 
in various astrophysical contexts. 



4.2 A framework for setting the reduced light 
speed value 

In the extremely complex framework of galaxy formation 
simulations, the accuracy of the results obtained using the 
RSLA can really only be assessed by convergence tests. It 
is nonetheless useful to consider a simple idealized setup in 
order to derive a physical intuition of where, when, and by 
how much one may reduce the speed of light. In this section. 



^ But ATON runs on GPUs, which are about a hundred times 
faster than CPUs, whereas RAMSES-RT runs on CPUs and thus 
can't afford such huge amount of RT subcycling. NB: ATON also 
increases the timestep by working on the coarse grid, and hence 
multiplying Ax by a factor ~ 2^-^ = 64 - 256 in Eq.lsTl 



we thus discuss the expansion of an ionized region around a 
central source embedded in a uniform neutral medium. 

We consider a source turning on and emitting ioniz- 
ing photons at a rate N into a homogeneous hydrogen-only 
medium of number density nn • An expanding sphere of ion- 
ized gas forms around the source and halts at the Stromgren 
radius rs within which the rate of recombinations equals the 
source luminosity: 



rs 



3iV 



V47ra^r 



1/3 



(32) 



where a^ ^ 
rate at T 



2.6 X 10~^^ cm^ s~^ is the case-B recombination 
^ 10^ K, and where we have assumed that the 
plasma within rs is fully ionized. 

The relativistic expansion of the I-front to its final ra- 



dius rs is derived in Shapiro et al. ( 2006| , and may be ex- 
pressed as: 



qy - ln(l 



y ), 



(33) 



where w — t/trec is time in units of the recombination time 
tree = (^hQ^^)~^, y = ri/rs is the position of the I-front 
in units of rs, and the factor q = tcross/trec = Ts/(ctrec) de- 
scribes the light crossing time tcross across the Stromgren ra- 
dius in units of the recombination time, and basically encom- 
passes all the free parameters in the setup (source luminos- 
ity, gas density, and temperature). Writing q oc N^'^n^ , 
we see that in many astrophysical contexts, q stays in the 
range ^ 10~^ — 10~^ (see Table fl]), simply because we are 
generally either interested in the effect of bright sources (e.g. 
a whole galaxy) on relatively low-density gas (e.g. the IGM) 
or of fainter sources (e.g. an 0-star) on high-density gas (at 
e.g. molecular-cloud densities). 

Let us now discuss briefly the evolution of an I-front 
given by Eq. |33]for three illustrative values of q\ 

• q — (blue curve of Fig. [^: This is the limiting non- 
relativist ic case, which assumes an infinite speed of light 
(tcross = 0). In this case, the I-front expands roughly as 



y (X. w ' (its speed decreases as w 



-2/3X 



almost all the way 



to rs, which it reaches after about a recombination time. 

• ^ = 1 (red curve of Fig. ^\ Here (and for all ^ > 1), the 
I-front basically expands at the speed of light all the way to 
rs, which it thus reaches after a crossing time. 

• q — 10~^ (green curve of Fig. IM): In this typical case, 
the I-front starts expanding at the speed of light, until w ~ 
(^/3)^/^. It then slows down and quickly reaches the limiting 
^ = behavior after a crossing time (at w ^ q). The I-front 
then reaches rs after a recombination time (at w ^ 1). 

An important feature appearing in the two latter cases is 
that for any physical setup q > 0, the F front is always 
well described by the q = limit after a crossing time (i.e. 
'^ ^ q)- We can use this feature to understand the impact of 
reducing the speed of light in our code. Say we have a physi- 
cal setup described by a value ^o • Reducing the speed of light 
by a factor fc<l{cr = fcc) implies an increase by a fac- 
tor 1/ fc of the effective crossing time, and the effective q in 
our experiment becomes qo/ fc- The solution we obtain with 
Or will be accurate only after an effective crossing time, i.e. 
after w = qo/fc- Before that time, the reduced- light-speed 
solution will lag behind the real one. 

How much one may reduce the speed of light in a given 
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Table 1. Stromgren sphere properties for typical cosmological regimes, with the inferred minimum ahowed hght speed fractions. 



Regime 


riH [cm 3] 


N [s-i] 


rs [kpc] 


Across [Myr] 


tree [Myr] 


Q 




Tsim [Myr] 


Wsim 


Jc,min 


MWISM 


10-1 


2 10^0 


0.9 


3 10-3 


1.2 


2 10- 


-3 


1 


1 


3 10- 


-2 


MW cloud 


102 


2 10^8 


2 10-3 


6 10-6 


1 10-3 


5 10- 


-3 


0.1 


80 


6 10- 


-4 


Iliev tests 1,2,5 


10-3 


5 1048 


5.4 


2 10-2 


122.3 


1.4 10 


-4 


10 


8 10-2 


2 10- 


-2 


Iliev test 4 


10-4 


7 10^2 


600 


2 


1200 


2 10- 


-3 


0.05 


4 10-5 


1 






-4 -3 -2 

Figure 2. I-front expansion in a Stromgren sphere for a set of 
values of the dimensionless crossing time q. The blue curve shows 
the infinite-light-speed limit {q = 0). The green curve shows a 
typical case with q = 10-3, and the red curve shows the q = 1 
case, as discussed in the text. The thin grey curves show other 
values of g, spanning the range 10-^ — 10 in steps of one dex. 
The grey lines in the top-left corner of the plot show slopes corre- 
sponding to an expansion at the speed of light (dot-dashed line) 
or as {t/trecY^^ (dashed line). For any g > 0, the I-front radius 
is accurately described by the q = 0-limit after a crossing time. 



numerical experiment then depends on the boundary con- 
ditions of the problem and their associated timescales. Call 
Tsim the shortest relevant timescale of a simulation. For ex- 
ample, if one is interested in the effect of radiative feedback 
from massive stars onto the ISM, Tsim can be set to the life- 
time of these stars. If one is running a very short experiment 



(see Sec. 6.5), the duration of the simulation may determine 
Tsim- Given this timescale contraint Tsim, one may reduce 
the speed of light by a factor such that the I-fronts will be 
correctly described after a timelapse well shorter than Tsim, 
i.e. tcross/fc ^ Tsim- In othcr words, one may typically use 
fc = min(l; ^ 10 x tcross/Tsim)- We now turn to a couple of 
concrete examples. 



4.3 Example speed of light calculations 

In Table IT] we take some concrete (and of course very ap- 
proximate) examples to see generally what values of fc are 
feasible. We consider cosmological applications from inter- 
galactic to inter-stellar scales and setups from some of the 
RT code tests described in ^ 



Reionization of the inter-galactic medium 

Here we are concerned with the expansion of ionization 
fronts away from galaxies and into the IGM, as for example 
in the fourth test of Iliev et al. (2006). In this test, the IGM 
gas density is typically nu = 10~^ cm~^, and the sources 
have iV = 7 10^^ s~^. In such a configuration, the Stromgren 
radius is rs '^ 600 kpc, corresponding to a crossing time 
Across ^ 2 Myr. Because of the low density of the gas, the 
recombination time is very long (> 1 Gyr), and we are thus 
close to the q = 10 ~^ case discussed above (the green curve 
inFig.[2|. 

Test 4 of 1106 is looked at output times Tsim,i =0.05 
Myr and Tsim, 2 = 0.4 Myr (see Fig. 19). In both cases, 
Tsim < tcross, and wc cauuot reduce the speed of light to get 
an accurate solution at these times, because the expanding 
front has not yet reached the q = limit. Interestingly, we 
cannot increase the speed of light either, as is done in 1106 
with C^-Ray which assumes an infinite light-speed. From 
Fig. [2] it is clear that this approximation (the blue curve) 
will over-predict the radius of the front. We can use the anal- 
ysis above to note that had the results been compared at a 
later output time Tsim > 2 Myr, the infinite-light-speed ap- 
proximation would have provided accurate results. It is only 
ten times later, however, that reducing the speed of light by 
a factor ten would have provided accurate results. 

We conclude that propagating an I-front in the IGM at 
the proper speed requires to use a value of the speed of light 
close to the correct value. This is especially true in Test 4 
of the Iliev et al. (2006) paper (last row of Table [l}. This 
confirms that for cosmic reionization related studies, using 
the correct value for the speed of light is very important. 



Inter-stellar medium 

There is admittedly a lot of variety here, but as a rough es- 
timate, we can take typical densities to be nu ^ 10~^ cm~^ 
in the large-scale ISM and nu ^ 10^ cm""^ in star-forming 
clouds. In the stellar nurseries we consider single OB stars, 
releasing Nob ^ 2 x 10^^ photons per second, and in the 
large-scale ISM we consider groups of {^ 100) OB stars. 
The constraining timescale is on the order of the stellar cy- 
cle of OB stars (Tsim '^ 10 Myr), and less for the stellar 
nurseries. In these two cases, which are representative of the 
dense ISM inside galactic disks, we see in Table [l] that the 
allowed reduction factor for the speed of light is much larger 
{fc — 10~^ to 10~^). This is due to two effects acting to- 
gether: the gas density is higher, but the sources are fainter, 
since we are now resolving individual stellar clusters, and 
not an entire galaxy. Test 1, 2 and 5 of Iliev et al. (2006) 
are also representative of such a favorable regime to use the 
reduced speed of light approximation (second to last row 
in Table [T| . This rigorous analysis of the problem at hand 
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confirms that propagating I-front in galaxy formation simu- 
lation can be done reliably using our current approach, while 
cosmic reionization problems are better handled with GPU 
acceleration and the correct speed of light. 



Normal (non- smoothed) RT 



4.4 Smoothed RT 

A problem we had to face, while performing RAMSES-RT 
galaxy formation runs, as well as the various test cases pre- 
sented here, is that there is often a small number of cells, 
usually along I- fronts, or close to strong radiation sources, 
that execute a huge number of thermochemistry subcycles in 
a single RT time- step. This is in part fault of the operator- 
splitting approach used, where the RT equations have been 
partly decoupled. Specifically, the photon density updates 
happen in three steps in this approach (see Fig.p] top). The 
photon injection step always increases the number of pho- 
tons, usually by a relatively large amount, and the transport 
step does the same when it feeds photons into cells along 
these I-fronts. The thermochemistry step in the I-front cells 
has the exact opposite effect: the photon density decreases 
again via absorptions. If the photon-depletion time is shorter 
than the Courant time, we have a curious situation where 
the cell goes through an inefficient cycle during the thermo- 
chemistry subcycles: it starts neutral with a large abundance 
of photons (that have come in via the transport and/or 
photon injection steps). It first requires a number of sub- 
cycles to evolve to a (partly) ionized state, during which the 
photon density is gradually decreased. It can then reach a 
turnaround when the photons are depleted. If the RT time- 
step is not yet finished, the cell then goes into a reverse 
process, where it becomes neutral again. This whole cycle 
may take a large number of thermochemical steps, yet the 
cell gas ends up being in much the same state as it started. 

In reality, the ionization state and photon density would 
not cycle like this but would rather settle into a semi- 
equilibrium where the rate of ionizations equals that of re- 
combinations. 

For the purpose of saving up on computing time and re- 
ducing the number of thermochemistry subcycles, we have 
implemented an optional strategy we call smoothed RT that 
roughly corrects this non-equilibrium effect of operator split- 
ting (see Fig. ^ bottom). In it, the result of (A/'/,F^) from 
the transport and injection steps in each cell is used to in- 
fer a rate for the thermochemistry step, rather than being 
set as an initial condition. We use the pre-transport, pre- 
injection values of Nj and F* as initial conditions for the 
thermochemistry, but instead update the thermochemistry 



equations (25) and (26) to 



dt 
dt 



Hl,Hel,Hell 



Y, UjCalN, + Nr + A^^, (34) 



y^ UjCo^jFi + Fi, 



(35) 



where the new terms at the far right represent the rates at 
which the photon densities and fluxes changed in the trans- 
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Many cooling steps! 









time 



Smoothed RT 



t + At 
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Figure 3. Sketch plots showing a photon density evolution over 
a global RT time-step with normal RT (top) and smoothed RT 
(bottom). In normal RT the photon density is updated to N' 
during photon transport (a) and injection (b). This is then used 
as an initial state for thermochemistry (c). It is often the case that 
the photons are depleted over the global time-step At, in a pro- 
cess which takes many thermochemistry subcycles. In smoothed 
RT, the photon density state is not updated by the transport and 
injection steps, but rather the difference is used to infer a pho- 
ton injection rate for the cell, which is gradually added during 
each thermochemistry substep. This can dramatically reduce the 
needed number of chemistry substeps. 



port and injection steps, i.e. 



N^ 



At 
F^- — F* 

^ I ■'- ^ 

At 



(36) 

(37) 



The idea is that when the photons are introduced like 
this into the thermochemistry step, they will be introduced 
gradually in line with the subcycling, and the photon density 
vs. ionization fraction cycle will disappear as a result and be 
replaced with a semi-equilibrium, which should reduce the 
number of subcycles and the computational load. The total 
photon injection (or depletion) will still equal Nl — Nf , so 
in the limit that there are no photoionizations or photon- 
emitting recombinations, the end result is exactly the same 
photon density (and flux) as would be left at the end of the 
transport and injection steps without smoothed RT. 

The gain in computational speed is dependent on the 
problem at hand, and also on the reduced light speed, which 
determines the size of the RT time-step. We've made a 
comparison on the computational speed between using the 
smoothed and non-smoothed RT in a cosmological zoom 



simulation from the NUT simulations suite (e.g. Powell, Slyz 
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fc Devriendt|20lT ) that includes the transfer of UV photons 
from steUar sources. Here, smoothed RT reduces the average 
number of thermochemistry subcycles by a factor of 6 and 
the computing time by a factor 3.5. So a lot may indeed be 
gained by using smoothed RT. 

One could argue that the ionization states in I-fronts 
are better modelled with smoothed RT, since the cycle of 
photon density and ionization fraction is a purely numeri- 
cal effect of operator splitting. We have intentionally drawn 
a slightly higher end value of Ni in the smoothed RT than 
non-smoothed in Fig.ls] whereas non-smoothed RT can com- 
pletely deplete the photons in a cell, smoothed RT usually 
leaves a small reservoir after the thermochemistry, that more 
accurately represents the "semi-equilibrium value" . 

Of course an alternative to smoothed RT, and a more 
correct solution, is to attack the root of the problem and 
reduce the global time-step length, i.e. also limit the trans- 
port and injection steps to the 10% rule. Reducing the global 
time-step length is highly impractical though; the main rea- 
son for using operator splitting in the first place is that it 
enables us to separate the timescales for the different steps. 

The same method of smoothing out discreteness that 
comes with operator splitting (in the case of pure hydrody- 
namics) has previously been described by Sun ( 1996), where 
it is referred to as "pseudo- non-time-splitting" . 



RADIATION HYDRODYNAMICS IN 
RAMSES 



RAMSES ( |Teyssier]|2002 ) is a cosmological adaptive mesh re- 
finement (AMR) code that can simulate the evolution and 
interaction of dark matter, stellar populations and baryonic 
gas via gravity, hydrodynamics and radiative cooling. It can 
run on parallel computers using the message passing inter- 
face (MPI) standard, and is optimized to run very large nu- 
merical experiments. It is used for cosmological simulations 
in the framework of the expanding Universe, and also smaller 
scale simulations of more isolated phenomena, such as the 
formation and evolution of galaxies, clusters, and stars. Dark 
matter and stars are modelled as collisionless particles that 
move around the simulation box and interact via gravity. 
We will focus here on the hydrodynamics of RAMSES though, 
which is where the RT couples to everything else. 

RAMSES employs a second-order Godunov solver on the 
Euler equations of gravito hydrodynamics in their conserva- 
tive form, 



dp 
di 



+ V • (pu) = 



dt 



(pu) + V • (pu u) + VP : 



dE 

dt 



+ V-((^ + P)u) 



(38) 
(39) 
pu- V(/) + A(p,£), (40) 



where t is time, p the gas density, u the bulk velocity, the 
gravitational potential, E the gas total energy density, P the 
pressure, and A represents radiative cooling and heating via 
thermochemistry terms (resp. negative and positive), which 
are functions of the gas density, temperature and ionization 
state. In RAMSES, collisional ionization equilibrium (CIE) is 
traditionally assumed, which allows the ionization states to 
be calculated as surjective functions of the temperature and 
density and thus they don't need to be explicitly tracked 




Figure 4. An oct - the basic grid element in RAMSES. 



in the code. E is divided into kinetic and thermal energy 
density (e) components: 



1 2 

-pu +£. 



(41) 



The system of Euler equations is closed with an equation of 
state which relates the pressure and thermal energy. 



P = (7 - l)e, 



(42) 



where 7 is the ratio of specific heats. The Euler equations 
are adapted to super comoving coordinates, to account for 
cosmological expansion, by a simple transformation of vari- 



ables (see J 5.4). 

The Euler equations are solved across an AMR grid 
structure. Operator splitting is employed for the thermo- 
chemistry source terms, i.e. A is separated from the rest 
of the Euler equations in the numerical implementation - 
which makes it trivial to modify the thermochemistry solver, 
i.e. change it from equilibrium to non-equilibrium. 

The basic grid element in RAMSES is an oct (Fig. l4|, 
which is a grid composed of eight cubical cells. A conserva- 
tive state vector U = (p, pu, E, pZ) is associated with each 
cell storing its hydrodynamical properties of gas density p, 
momentum density pu, total energy density E and metal 
mass density pZ. (One can also use the primitive state vec- 
tor, defined as W = (p, u, P, Z).) Each cell in the oct can 
be recursively refined to contain sub-octs, up to a maximum 
level i of refinement. The whole RAMSES simulation box is 
one oct at ^ = 1, which is homogeneously and recursively 
refined to a minimum refinement level ^min, such that the 
coarse (minimum) box resolution is 2^"^^" cells on each side. 
Octs at or above level imin are then adaptively refined during 
the simulation run, to follow the formation and evolution of 
structures, up to a maximum refinement level ^max, giving 
the box a maximum effective resolution of 2^'^'^'^ cell widths 
per box width. The cell refinement is gradual: the resolu- 
tion must never change by more than one level across cell 
boundaries. 



5.1 RAMSES multi-stepping approach 

With AMR multi- stepping, the resolution is not only adap- 
tive in terms of volume, but also in time, with different 
timestep sizes on different refinement levels. A coarse time- 
step, over the whole AMR grid, is initiated at the coarse 
level, ^min, as we show schematically in Fig. |5] First, the 
coarse time-step length Ati^^^ is estimated via (the min- 
imum of) Courant conditions in all ^min cells. Before the 
coarse step is executed, the next finer level, imin + 1, is 
made to execute the same time- step, in two substeps since 
the finer level Courant condition should approximately halve 
the time- step length. This process is recursive: the next finer 
level makes its own time-step estimate (Courant condition. 
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Figure 5. Recursive hydro time-stepping over one coarse time- 
step in the AMR levels of RAMSES, here shown for a three-level 
AMR structure. Each solid arrowed line represents a time-step 
which is executed for all cells belonging to the corresponding 
AMR level. The numbers indicate the order of the time-stepping, 
including the calls to finer levels (1, 2, and 6). 







Listing 1: The 


AME 


I step in RAMSES. 




recursive subroutine 


amr. 


.step(^): 






if t 


< 


^max and any ce 


lis 


exist in 


^+1 




ca 


11 


amr_step (-^ + 1) 










ca 


11 


amr_step (^ + 1) 










call 


h 


ydro_solver {i) : 


all 


i cells 


and some ^ — 1 


call 


e 


q_thermochemist 


ry(^): all i 


leaf 


cells 


end 















but also At^ ^ At^-i) and has its next finer level to ex- 
ecute two substeps. This recursive call up the level hier- 
archy continues to the highest available level ^max, which 
contains only leaf cells and no sub-octs. Here the first two 
substeps are finally executed, with step lengths At^^^^ ^ 
^^^min/2^"'''''"^"''''- When the two -£max substeps are done, 
the ^max — 1 time-step is re-evaluated to be no longer than 
the sum of the two substeps just executed at ^max, and then 
one -£max — 1 step is executed. Then back to level -£max to 
execute two steps, and so on. The substepping continues in 
this fashion across the level hierarchy, ending with one time- 
step for the coarsest level cells (with a modified time-step 
length At,^^ J. 

At the heart of RAMSES lies a recursive routine called 
amr_step(-£) which describes a single time-step at level ^, 
and is initially called from the coarsest level (^min)- To fa- 
cilitate our descriptions on how the RT implementation is 
placed into RAMSES, we illustrate the routine in pseudocode 
format in Listing [l] where we have excluded details and bits 
not directly relevant to RHD (e.g. MPI syncing and load- 
balancing, adaptive refinement and de-refinement, particle 
propagation, gravity solver, star formation, and stellar feed- 
back) . 

First, the recursion is made twice, solving the hydrody- 
namics over two sub-steps at all finer levels. Then the Euler 
equations are solved over the current coarse time-step, for all 
cells belonging to the current level. It is important to note 
here that the hydrodynamical quantities are fully updated 
at the current level in the hydro_solver, but there are also 



intermediate hydro updates in all neighboring cells at the 
next coarser level. The coarser level update is only partial 
though, because it only reflects the intercell fluxes across 
inter-level boundaries, and fluxes across other boundaries 
(same level or next coarser level) will only be accounted 
for when the coarser level time- step is fully advanced. Until 
then, these coarser level neighbor cells have gas states that 
are not well deflned, since they only reflect some of their in- 
tercell fluxes. It effectively means that at any point between 
the start and flnish of the primary (coarse) call to ainr_step, 
there are some cells in the simulation box (lying next to flner 
level cells) that have ill-deflned intermediate hydrodynami- 
cal states. This point is further illustrated in Appendix [D] It 
is important to keep in mind when considering the coupling 
of RT with the hydrodynamics of RAMSES. 

Having put down the basics of AMR hydrodynamics, 
we are now in a position to add radiative transfer. 

5.2 RAMSES-RT 

In RAMSES-RT, each cell stores some additional state vari- 
ables. Here W = (p, pU, £;, yoZ, pXHll, pXHell, P^Helll, A^i, F^), 

where xhii, x^eii and XHeiii are the hydrogen and helium ion- 
ization fractions, which are advected with the gas as passive 
scalars (in the hydro solver), and A/'^, F^ represent the AM 
variables of photon density and flux for each of the M pho- 
ton groups. Note that this represents a hefty increase in the 
memory requirement compared to the hydrodynamics only 
of RAMSES: the memory requirement for storing lA (which 
is the bulk of the total memory in most simulations) is in- 
creased by a factor of 1.5(l+4/9M), where the 1.5 represents 
the ionization fractions and the parenthesis term represents 
the photon fluxes and densities. Thus, with three photon 
groups, the memory requirement is increased by roughly a 
factor 3.5 compared to a traditional RAMSES simulation. 

Given the time-scale difference between hydrodynamics 
and radiative transfer, the obvious approach to performing 
RHD is to sub-cycle the three radiative transfer steps (in- 
jection, advection, thermochemistry) within the hydrody- 
namical step. There is, however, a major drawback to this 
approach, which is that it is incompatible with AMR multi- 
stepping: the RT sub-cycling must be done before/after each 
hydrodynamical AMR step at the finest refinement level 
only, and since light can in principle cross the whole box 
within the flne level hydrodynamical timestep, the RT sub- 
cycling must be done over the whole grid, over all levels. 
However, the partial hydrodynamical flux between cells at 
level boundaries always leaves some cells between the fine 
level steps with an intermediate (i.e. partially updated) gas 
state. This makes the thermochemistry ill-defined in those 
cells, since it needs to update the gas temperature in every 
cell, and for this to work the temperature must have a well 
defined and unique value everywhere. There are three ways 
around this: 

First is to perform the RT subcycling only after a coarse 
hydrodynamical step, but here potentially thousands of fine- 
scale hydro steps would be executed without taking into 
account the thermochemistry. 

Second, is to prohibit AMR multi-stepping, which 
makes the whole grid well defined after each step and thus 
allows for RT sub-cycling over the whole box. Multi-stepping 
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Figure 6. Diagram of the ainr_step in RAMSES-RT. This is much 
hke the normal amr_step in RAMSES, except that the time-step 
length has the extra constraint of the light speed Courant con- 
dition, and each level i step also performs photon injection, RT 
transport and thermochemistry over the same time-step and level. 



is however one of the main advantages of AMR, and essen- 
tially allows us to refine in time as well as space, so this isn't 
really an option. 

We thus default to the third strategy, which we use 
in RAMSES-RT. Here we drop the subcycling of RT within 
the hydro step and perform the two on the same timestep 
length, which is the minimum of the RT and hydro timestep. 
Thus, with each hydro step, at any level, the RT steps are 
performed over the same level only. The basic scheme is 
illustrated in Fig. [6] and the pseudocode for the updated 
cunr.step is shown in Listing l2] Obviously, the main draw- 
back here is the timescale difference, which can be something 
like a factor of 100 — 1000, meaning the number of hydrody- 
namical steps is increased the same factor and the run-time 
accordingly (plus numerical diffusion likely becomes a prob- 
lem with such small hydrodynamical time-steps). However, 
if we also apply a reduced speed of light, we can shrink 
this factor arbitrarily, down to the limit where the hydro- 
timestep is the limiting factor and the only increase in com- 
putational load is the added advection of photons (which is 
considerably cheaper for one photon group than the hydro- 
dynamical solver) and the non-equilibrium thermochemistry 
(which typically has a computational cost comparable to the 
equilibrium solver of RAMSES, provided we use RT smooth- 



ing). The question, which we have tried to answer in [4.2 



is then how far we are allowed to go in reducing the light 
speed. 

Parallelization is naturally acquired in RAMSES-RT by 
simply taking advantage of the MPI strategies already in 
place in RAMSES. 



5.3 Radiation transport on an AMR grid 

In RAMSES-RT, the radiation variables are fully incorporated 
into the AMR structure of RAMSES. The ionization fractions 
and photon densities and fluxes are refined and de-refined 
along with the usual hydro quantities, with a choice of inter- 
polation schemes for newly refined cells (straight injection 
or linear interpolation). The radiative transfer, i.e. injection, 
transport and thermochemistry, is multi-stepped across the 
level hierarchy, thus giving AMR refinement both in space 
and time. Inter-level radiation transport is tackled in the 





Listing 2: The AMR step in RAMSES-RT. 


recursive subroutine amr_step(^) 


if i 


< -^max and any cells in i -\- 1 


ca 


11 amr_step(-^+ 1) 


ca 


11 amr_step(-^+ 1) 


call 


photon_injection_step (i) 


call 


hydro-solver (^) : all i cells and some i—1 


call 


rt -transport (-^) : all i cells and some i—1 


call 


neq_thermochemistry (-^) : all £ leaf cells 


end 





same way as the hydrodynamical advection, i.e. transport 
on level i includes partial updates of neighbouring cells on 
level i—1. Update of the finer level cell RT variables over 
level boundaries involves the RT variables in a coarser cell, 
which are evaluated, again with the same choice of inter- 
polation schemes. RAMSES-RT includes optional refinement 
criteria on photon densities, ion abundances and gradients 
in those, in addition to the usual refinement criteria that 
can be used in RAMSES (on mass and gradients in the hydro- 
dynamical quantities). 

Of the standard RT and RHD tests described in 
Section [6] four include active or inactive grid refine- 
ment, demonstrating that the radiation hydrodynamics 
perform robustly in conjunction with (on-the-fly) cell 
refinements/de-refinements. In addition, we demonstrate in 
Fig. |7| how radiation flux is well retained across changes in 
grid refinement. The upper left map of the figure shows a 
beam of radiation in a 2D RAMSES-RT experiment, where 
we use the HLL flux function and deactivate radiation-gas 
interactions. The beam is injected into two cells in the bot- 
tom left corner by imposing a unity reduced photon flux of 
3 X 10^*^ photons s~^, corresponding to a photon density of 
1 cm~^, at an angle of 26.5° from the horizontal. The beam 
traverses a circular region of 2 successive levels of increasing 
refinement, going from refinement level 6 to 8, i.e. effective 
resolutions of 64^ to 256^ cells. We use here straight injection 
(i.e. no interpolation) for inter-level cell fluxes, but linear in- 
terpolation gives identical results. To the right of the map 
we plot photon flux profiles, cN, across the coloured ver- 
tical lines in the map. The beam experiences diffusion, but 
this is exclusively due to the intercell flux function and inde- 
pendent of the refinement changes. The far left plot shows 
the integrals across each flux profile, i.e. the total photon 
flux across each line. The values are consistent to about 1 
in lO'^, and actually the only visible discrepancy appears 
to be a boundary effect close to the right side of the box 
(the boundaries allow free outflow). The total flux is thus 
clearly very well preserved across the grid refinement levels. 
To further demonstrate this, the lower panel in the same 
figure shows an identical experiment except that the beam 
is horizontal, such that it can be perfectly maintained with 
the HLL flux function. There is no diffusion whatsoever re- 
sulting from the change in refinement along the path of the 
beam, as can be seen from the map and from the flux profiles 
across the vertical lines in the plots (the profiles are plot- 
ted with increasing thickness to show that they coincide). 
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Figure 7. 2D beam experiments, demonstrating photon flux conservation across changing reflnement levels. The upper panel shows an 
experiment with an off-axis beam, 26.5° from the horizontal, and the lower panel is an identical experiment except the beam is horizontal. 
The maps on the right show photon number density, with the grid structure overlaid in grey (which is kept constant throughout the 
experiments). Coloured vertical lines mark x-positions at which photon flux proflles are plotted in the left plots. The right plots show 
integrals of each proflle, i.e. the total photon flux across each x-coordinate. 
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Figure 8. 2D beam experiment, same as Fig. [7| but with on-the-fly AMR reflnement. 



The total flux is retained perfectly to the number precision, 
which here is 7 decimals. 

We also consider another beam with the same setup, 
shown in Fig. [S] where instead of a static refinement region, 
the grid is actively refined on inter-cell gradients in photon 



density N. According to the criterion, two adjacent cells at 
positions i and i -\- 1 are refined if 



>0.4. 



(43) 



m + AT^+i + 10-3 cm-2 
Straight injection (no interpolation) is used here for inter- 
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level fluxes and cell refinements, but the results are identical 
when linear interpolation is used for inter-level fluxes and 
cell refinements. The snapshot here is taken shortly before 
a light crossing time, in order to check for flux conservation 
close to the far end of the beam where the grid is being 
actively refined. The plot on the far right shows the flux 
conservation across different x-coordinates. (Note the total 
flux is slightly different from that in Fig. [7| because of the 
different geometry of the beam injection.) The total flux is 
maintained to the order of 0.25%, except for the outermost 
profile (red cross), where the relatively large discrepancy is 
due to the placement close to the end of the beam. If the 
experiment is let to run for more than a light crossing time, 
the discrepancy levels out to within 0.03% of the innermost 
(black cross) profile. 

These simple beam experiments demonstrate that the 
code accurately transports radiation across (even dynami- 
cally) changing refinement levels. The main errors are the ar- 
tificial diffusion of radiation on the grid, which is not caused 
by refinement, but rather by the inter-cell flux function, and 
the dipole approximation inherent to the Ml closure, which 
does not allow opposing streams of radiation to pass through 
one another. Note though that while the diffusion is artifi- 
cial, the total flux is well maintained, i.e. energy is conserved. 



5.4 Cosmological settings 

RAMSES uses super-comoving variables to allow for the im- 
pact of the cosmological expansion on the Poisson equation, 
the equations of hydrodynamics ( 38][4Q ) and particle prop- 
agation (Martel & Shapiro 1998| Teyssier 2002): a change 



is made from the physical variables to super-comoving ones 
with 



dt = -^r- dt, 



aL 



HoL 



^^mPcJ^Q J^ 



1 LmPc-tlQ 



where Hq is the Hubble constant, Qm the matter density pa- 
rameter, L the comoving width of the simulation box (phys- 
ical width at a = 1), and pc the critical density of the Uni- 
verse. When these variables are used instead of the physical 
ones, the cosmological expansion is accounted for, while all 
relevant equations remain unchanged, Euler equations in- 
cluded. 

For consistency, and to partly account for the effect of 
cosmological expansion on the radiative transfer, the addi- 
tional change is made in RAMSES-RT to super-comoving RT 
variables for the photon transport: 



N = a^ N, 



F = 



HoL 



HoL 



The dilution (oc a~^) of photon number density is thus ac- 
counted for, while it can easily be verified that Eqs. (|6|-([7|) 
remain unchanged with the new variables - including the 
reduced flux (14) used in the Ml tensor (12). 



Note that when reduced light speed is used, the pho- 
tons will be over-diluted in cosmological simulations, since 
the time taken for them to get from source to destination 
will be overestimated. Note also that wavelength stretching 
with redshift, which in reality adds a fourth power of a to 



the dilution of N^, is not accounted for here. This is actu- 
ally non-trivial to do: one could add one power of a to the 
definitions of N and F, but it would be a very crude ap- 
proximation of the wavelength dilution, as the wavelength 
shift that should feed photons from one group to the next is 
neglected. In any case, this effect is likely to be important 
only in the context of reionization, where the photons have 
a chance of travelling cosmological distances before they are 
absorbed. To the best of our knowledge, at this time no pub- 
lished RT or RHD simulations account at all for the effect 
of cosmological expansion on the photon transport. 



6 RADIATIVE TRANSFER TESTS 

The tests described in this section come from two papers 
that were born out of a series of workshops on radiative 
transfer. Tests with simple analytic results to compare to 
are hard to engineer in radiative transfer, so the solution 
was to instead make simple tests where the correct result is 
not necessarily well known but the results of many different 
codes can instead be compared. Thus it is likeliest that the 
correct results are usually where most of the codes agree, 
and if a code stands out from all or most of the others in 
some way, this would most likely be a problem with that 
particular code. These tests have become sort of benchmark 
tests for RT codes, and most publications that present new 
implementations use some or all of these tests for validation. 



The first paper is Iliev et al. (2006a), hereafter known 



as 1106 - it describes four RT post-processing tests, i.e. with 
the hydrodynamic advection turned off, and shows the re- 



sults for 11 RT codes. The second paper is Iliev et al. (2009 ), 
hereafter known as 1109 - it describes three additional tests, 
and results for 9 codes, where the RT is coupled to the hy- 
drodynamics. 

The tests results from 1106 and 1109 are normally down- 
loadable on the web, but at the time of this writing the links 
have been down for some time. However, Ilian Iliev has been 
kind enough to provide all test results for one of the codes, 
the grid based short characteristics code C^-Ray, which is de- 
scribed in detail in Mellema et al. (2006). We thus present 



here RAMSES-RT results with comparisons to those of C -Ray. 
The inclusion of the C^-Ray results in the plots shown here 
should be useful to guide the eye if one then wants to com- 
pare with the other codes in 1106 and 1109. 

As prescribed by the test papers, all tests use hydrogen 
only gas. We use smooth RT in the RAMSES-RT runs for all 
tests, but remark that turning off the smoothing has no 
discernible effect on the results (only calculation speed). 



6.1 



1106 test 0: 
physics 



The basic thermochemistry 



This is essentially a one-cell test of the non-equilibrium ther- 
mochemistry and not radiative transfer per se, so it doesn't 
really count with the rest of the comparison project tests 
(hence test zero). It is important nontheless since thermo- 
chemistry is a major new component in RAMSES-RT. 

We start with completely neutral hydrogen gas with 
density uh — 1 cm~^ and temperature T — 100 K at t = 0. 
A photo-ionizing flux of F = 10^^ s~^ cm~^ with a 10^ K 
blackbody spectrum is applied to the gas and maintained 
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Figure 9. 1106 test 0. Single-zone photoheating and ionization 
with subsequent cooling and recombinations. 




Figure 10. 1106 test 1. Map of the neutral fraction in a box slice 
at 2; = 0, at 500 Myr. Overplotted is the AMR grid, which is 
refined on the fly during the experiment from 64"^ to 128^ cells ef- 
fective resolution. Maximum reflnement stays on the corner source 
throughout the run, and it adaptively follows the expansion of the 
I-front. 



until t = 0.5 Myr at which point it is switched off. The run 
is continued for a further 5 Myr, allowing the gas to cool 
down and recombine. The resulting evolution of the neutral 
fraction and temperature of the gas is shown in Fig. [9] The 
evolution closely follows that of the codes described in 1106, 
with the exception of Simplex and FFTE which stand out 
somewhat, and we don't see any sign of the stiffness-induced 
oscillations that can be seen in the Crash code test. 



6.2 1106 test 1: Pure hydrogen isothermal HII 
region expansion 

A steady monochromatic (hv — 13.6 eV) source of radiation 
is turned on in a homogeneous neutral gas medium, and 
we follow the resulting expansion of a so-called Stromgren 
sphere of ionized gas. Heating and cooling is turned off and 
the temperature is set to stay fixed at T = 10^ K. The OTSA 
is applied in this test and we use GLF intercell fluxes. 

The box is a cube of width Lhox = 6.6 kpc. The gas 
density is nn = 10"*^ cm~^ and the initial ionization frac- 
tion is xhi = 1.2 X 10""^, corresponding to collisional ion- 
ization equilibrium. The radiative source is in the corner of 
the box and the emission rate is iV-^ = 5 x 10^^ photons s~^. 
We use the RSLA in this test with the light speed fraction 
set to /c = 10~^. AT08 contains an analysis of the effect 
of different light speeds in this test and others from 1106, 
and finds the results start diverging non-negligibly some- 
where between /c = 10~^ and 10~^. The simulation time 
is tsim = 500 Myr. As with all the tests in 1106 and 1109 
(except for test which doesn't involve a volume), the pre- 
scribed grid resolution is 128^ cells, but here we do things 
slightly differently to demonstrate on-the-fly AMR at work. 
Our box resolution is 64"^ cells, but we allow for one level 
of further refinement, i.e. to an effective resolution of 128^ 
cells that matches the prescribed resolution. The refinement 
criterion is applied on gradients in xhi and xhii- According 
to it, two adjacent cells at positions i and i + 1 are refined if 



X* + x^+i 



where x is either xhi or xhh- 



>0.8, 



(44) 



The Stromgren radius, rs, is the radius of the ionization 
front (I-front) from the center when steady state has been 
reached, and in the case of fixed density and temperature 
it has the simple analytical result shown in Eq. |32[ In this 
result the I-front evolves in time according to 



r\ 



rs 



[1 - e-'^'-] 



1/3 



(45) 



where tree = (^h^hii)""^ is the recombination time. For the 
parameters of this experiment, tree = 122.4 Myr and rs = 
5.4 kpc. 

Fig. [To] shows maps at 500 Myr of the neutral fraction, 
with the grid refinement overplotted, in a box slice at z = 0. 
The Stromgren sphere is nicely symmetric and qualitatively 
it can be seen to agree well with results from the RT codes 
described in 1106 (their Fig. 6). 

Fig. |lla| shows the evolution of the I-front position and 
velocity with RAMSES-RT (solid blue), compared with the an- 
alytic expression (green dot-dashed) and the result for the 
C^-Ray code (red dashed), which is typical for the RT code 
results presented in 1106 and does not stand out particularly 
in this test. Our result can be seen to match the C^-Ray one, 
though we have an initial lag due to the reduced speed of 
light that can best be seen in the top plot showing the frac- 
tion of the numerical result's I-front radius versus rs- The 
analytic ri is typically ahead of rs by < 5%, which is simply 
because the analytic result is step-like with complete ioniza- 
tion within rs and none outside, whereas the real result has 
a gradually evolving ionization profile with radius. Indeed, 
Pawlik & Schaye ( 2008| computed the exact analytic result 
to this problem, accounting for an equilibrium neutral frac- 
tion inside the Stromgren sphere, and found an equilibrium 
I-front radius which is exactly 1.05 rs. 

Fig. |llb| shows spherically averaged radial profiles of the 
gas ionization state at 30 and 500 Myr. Again we see a good 
match with the C^-Ray result. There is still a little lag in 
the I-front position at 30 Myr due to the RSLA and xhi is 
somewhat lower inside the Stromgren sphere in RAMSES-RT. 
However, the C^-Ray result stands out a little in this test in 
1106 as being most effective at ionizing the gas within the 
Stromgren sphere (i.e. has the lowest values of xhi), and the 
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Figure 11. 1106 test 1. (a) Evolution of the I-front position and velocity. Blue solid lines show our result, red dashed lines show the 
C^-Ray result and green dot-dashed lines show the analytic expression, (b) Spherically averaged profiles for neutral fractions xui and 
ionized fractions xhii at 30 and 500 Myr versus radius (in units of the box width Lbox)- (c) Histogram showing fractions of cells within 
bins of xiii at three simulation times. 
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Figure 12. 1106 test 1. Evolution of the globally averaged neutral 
fraction. 



RAMSES-RT result is firmly typical of the 1106 codes' results 
in this plot. 

A further comparison is made in Fig. |llc| here compar- 
ing ionization fraction histograms at three simulation times. 
Again the RAMSES-RT result closely matches the C^-Ray one, 
whose histograms fall into a group with the codes IFT, 
Flash-HC and FFTE that stand out a little in 1106 (Fig. 9) as 
having less frequent intermediate neutral fractions than the 
other codes. 

Finally for this test. Fig. [12] shows a comparison with 
C^-Ray of the globally averaged neutral fraction as a function 
of time. It is a close match, and the C^-Ray result is here 
firmly typical for the 1106 codes. 

All in all, there is nothing out of the ordinary in the 
RAMSES-RT result for 1106 test 1, except for a slight initial 
delay of the I-front which is to be expected due to the RSLA. 



6.3 1106 test 2: HII region expansion and the 
temperature state 

The setup here is the same as in 1106 test 1, except for the 
following points: 



• We allow for cooling and photo-heating of the gas, i.e. 
the temperature is no longer constant, and the analytic re- 
sult, Eq. 32 no longer applies (because of the non-constant 



recombination rate). 

• The initial temperature is 100 K. 

• The initial ionization fraction of the gas is xhh = 10~^. 
It should be fully neutral according to the test recipe, but 
this is (the default) minimum value for xhii in RAMSES-RT, 
that exists in order to keep bounds on the subcycling of the 
thermochemistry. 

• The radiation source is modelled as a blackbody with 
temperature T = 10^ K. The emission rate is the same as 
before, iV^ = 5 x 10^^ photons s~^ 

• Multi-frequency is approximated with three photon 
groups of Hi, Hei and Hell ionizing photons, as in all fol- 
lowing tests. The light speed fraction is fc = 10~^. 

• We don't use grid refinement in this test. The grid is 
homogeneous and the resolution is 128^ grid cells, as pre- 
scribed in 1106. 

Slice maps at z = of the neutral fraction and temper- 
ature are shown in Fig. |13| Both the ionization and heating 
fronts are smooth and symmetric, and the maps agree quali- 
tatively with other codes in 1106 (Figs. 11-14). In comparison 
with the same test with ATON (AT08, Fig. 3), both fronts 
are clearly much thicker here, which is due to our multi- 
frequency implementation (whereas ATON used one photon 
group). More detailed comparison with the 1106 codes can 
be ma de th rough the ionization state and temperature plots 
where we include the C^-Ray result. The ioniza- 



14a 



in Fig. 

tion state profile develops very similarly to that of C^-Ray, 
though we have less ionization on both sides of the front, 
especially on the outer side where the difference in xhii is 
as high as a factor of ten. Presumably this is due to the 
different implementations of multi-frequency photo-heating 
and cooling. The thermal profiles are also similar to C^-Ray, 
though we have considerably lower (up to a factor of two) 
temperatures on the inside of the I-front, and conversely 
higher temperatures on the outside. As can be seen in Fig. 
17 in 1106, C^-Ray has the strongest heating of any code on 
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Figure 13. 1106 test 2. Maps showing slices at 2; = of the neutral fraction and temperature at 10 Myr and 100 Myr. 
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Figure 14. 1106 test 2. (a) Evolution of the temperature and ionization state profiles, (b) Evolution of the ionization front. The top 
plot shows the ratio of the radius of the I-front in the tests, rnum versus the time-evolving radius ranalyt i^i the analytic result from test 
1 (Eg. |45| ). The middle plot shows the ratio of the test I-front radius versus the steady-state radius in the same analytic result (Eg. |32| ). 
The bottom plot shows the speed of the I-front, vi in units of a 'characteristic' speed, given by rs/^rec- (c) Histograms of temperature 
and ionized fraction. 



the inside of the I-front in this test and most codes have 
stronger heating on the outside, so our thermal profiles (as 
the ionization state profiles) are fairly typical of the ones 
presented in 1106 for this test. 

Fig. |14b| shows the evolution with time of the ioniza- 
tion front, compared with C^-Ray and the analytic result 
from test 1. The front moves more slowly here than in test 
1 due to the lower initial temperature, so we no longer lag 
behind in the initial front propagation. Our front propagates 
slightly further than in C^-Ray, and ends at almost exactly 
the same radius as the FFTE code, which has the furthest 
expanding I-front of any code in this test in 1106. Still the 
difference between the codes is small, with the ratio between 
the numerical and analytic results (rnum/Tanaiyt) ranging be- 
tween 1.01 and 1.11. 

Fig. |14c| shows histograms of the ionized fraction and 
temperature at different times in the test for RAMSES-RT 
and C^-Ray. The ionized fraction histograms are quite simi- 
lar, the biggest difference being a higher fraction of almost 
completely neutral gas xhh % 10~^ in RAMSES-RT, which we 
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Figure 15. 1106 test 2. Time evolution of the volume average 
neutral fraction. 



already saw in Fig. 14a (top) beyond the I-front. The tem- 
perature histogram for RAMSES-RT differs a bit from C^-Ray 
in having less extreme temperatures (C^-Ray has both hot- 
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ter gas and colder gas) but are very similar to those for the 
codes ART, RSPH and Crash in 1106. 

Finally, Fig. |15| shows the time evolution of the volume 
averaged neutral fraction in RAMSES-RT and C^-Ray, and here 
we see a close match. There is quite a lot of discrepancy 
between the different codes in the analogue plot in 1106 (Fig. 
20), with 3 groups of results, and our result closely follows 
those of C^-Ray, Crash and RSPH. 

As with test 1, there is nothing out of the ordinary in 
the RAMSES-RT result for 1106 test 2, except perhaps for an 
ever so slightly further advanced I-front than most codes in 
1106 have. 



6.4 1106 test 3: I-front trapping in a dense clump 
and the formation of a shadoAv 

This test considers self-shielding within a dense gas cloud 
bombarded on one side by UV radiation, and the shadow 
trailing on the 'dark' side - something which may find place 
with clouds close to sites of star- format ion. 

The setup is as follows: the simulation box has width 
Lho^ — 6.6 kpc. We place a spherical cloud of gas in the 
center of the (2/, ^)-plane, with radius rdoud = 0.8 kpc, and 



it's center at {xc^Vc^Zc) — (5,3.3,3.3), as seen in Fig. 16 
top left, showing an (x,^)— slice of gas density through the 
middle of the box. Outside the gas cloud we have n^^ 



T°^' = 8000 K and x^fi = 0, and inside 
= 200 ng"' = 4 X 10-^cm-^ T^^°^^ = 40 
10~^. We apply a constant ionizing photon 
^ cm~^ from the x — boundary of the box 



2 X 10"^ cm"*" 
we have n^^°"^ 
K and x^^°"^ = 
flux F = 10^ s" 
(left in the Fig. 16 maps), and run for 15 Myr. The spec- 
trum is a blackbody with effective temperature Teff = 10^ 
K, and as before we approximate mult i- frequency by split- 
ting into three photon groups bordered by the Hi, Hei, Hell 
photoionization frequencies. We use a light speed fraction of 
fc — 10~^. As usual the resolution prescribed by 1106 is 128^ 
cells, but here we apply static AMR refinement such that the 
coarse resolution is 64"^ cells, but a rectangular region that 
encompasses the gas cloud and the shadow behind it has 
one level of additional refinement, making the effective reso- 
lution in the cloud and its shadow 128^ cells. The refinement 
region is shown in the top panel of Fig. [16] plotted over a 
density map that shows the spherical gas cloud. In order to 
best capture the formation of a shadow behind the cloud, 
we apply the HLL flux function in this test rather than the 
usual GLF function, and we use the OTSA. We have run 
identical tests though, one with the GLF flux function, and 
one where we use the HLL flux function but don't assume 
the OTSA, and we show maps of those experiments for a 
qualitative comparison. 

Fig. [16] shows slices at z = 0.5 Lbox of the neutral frac- 
tion and temperature at 1 and 15 Myr. From second top to 
bottom row are shown RAMSES-RT+HLL, RAMSES-RT+HLL 
without the OTSA, RAMSES-RT+GLF (with the OTSA) and 
C^-Ray. The Lfront travels fast through the diffuse medium 
outside the cloud, but moves much more slowly inside it, and 
a shadow is cast behind it. As the UV radiation eats its way 
into the cloud, ionizing and heating it, the shadow also very 
slowly diminishes in width because some photons manage to 
cross through the edges of the cloud. The RAMSES-RT+HLL 
maps compare very well with C^-Ray, though the shadow 
is slightly thinner at 15 Myr and there is stronger heating 



inside the shadow; this could be due to differences in the 
multifrequency approach and/or photoheating. Without the 
OTSA, the shadow is diminished from the sides due to pho- 
tons being cast from the surrounding gas. Using the GLF 
flux function has much the same effect as not assuming the 
OTSA, though the shadow is considerably more diminished 
here. The result with HLL but without the OTSA is the 
most physical of the RAMSES-RT results, as one should ex- 
pect recombination photons to be cast into the shadow. 

Fig. |17a| shows the evolution of the position and speed 
of the Lfront through the center of the (2/, z)— plane. In solid 
blue we plot the RAMSES-RT result and in dashed red is the 
C^-Ray result for comparison. Horizontal dotted lines mark 
the edges of the cloud. There is a large initial delay in the 
Lfront compared to C^-Ray, which is because in the diffuse 
gas outside the cloud, the Lfront speed is limited by the 
reduced speed of light. After the Lfront gets into the cloud 
(lower dotted line) it quickly catches up and then evolves 
in a similar fashion in the two codes. If compared to the 
rest of the codes in 1106, it turns out that the evolution of 
the I-front in C^-Ray slightly stands out from the rest of 
the codes (e.g. a small upwards 'bump' in the front position 
at log {t/ tree) "^ 0, and a slightly shorter distance of the I- 
front from the origin at the end of the simulations), and 
most of the others in fact evolve very similarly to that of 
RAMSES-RT. The comparison appears best with RSPH, which 
has the furthest extended I-front at the end-time of 15 Myr. 
The same can be said for the speed of the front. If we look 
away from the initial ^ 0.2 Myr, when our I-front has to 
catch up, the speed compares reasonably to C^-Ray, and 
quite well to the other codes in 1106. 

Fig. 1 17b] shows the evolution of the mean ionized frac- 
tion and temperature inside the cloud, compared between 
RAMSES-RT and C^-Ray. The evolution is similar between the 
two codes in both cases. Compared with the other codes in 
1106, the evolution of the ionized fraction is most similar to 
RSPH, IFT and Coral, while the temperature in RAMSES-RT 
is consistently a little higher than in most codes (all except 
Coral and Flash which stand out quite a lot in mean tem- 
perature) . 

Fig. |17c| shows profiles of the ionization state and tem- 
perature along the x-axis at the center of the {y, z)— plane at 
1, 3 and 15 Myr. The ionization state profile in RAMSES-RT is 
similar in most respects to that of C^-Ray, though it extends 
a bit further at the end of the run-time. There is initially 
less ionization on the far side of the front in RAMSES-RT, 
but at the end of the run this is reversed and we have 
slightly more ionization on the far side in RAMSES-RT. This 
'shift' can be explained by the temperature profiles: at early 
times the cloud is efficiently shielding the far side from even 
the high-energy photons in both codes, but at the end of 
the RAMSES-RT run the shielding buffer in the cloud is thin 
enough that the high-energy photons can get through, hence 
efficiently heating the gas inside the buffer as well as in the 
shadow, and the gas in the shadow becomes slightly ionized 
as a consequence. The analogue ionization state profiles for 
the other codes in 1106 are mostly similar to ours. Most of 
them are actually closer to the RAMSES-RT than the C^-Ray 
profile, with the exception of Crash which has a much more 
underdeveloped I-front and less ionization, and FFTE and IFT 
which have an almost step- wise xnii-profile on the far side of 
the I-front. The temperature profiles differ pretty widely be- 
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Figure 16. 1106 test 3. Maps showing slices dX z = 0.5 Lbox- The top map shows the (constant) density field, with the static refine- 
ment overplotted. The second row shows the RAMSES-RT+HLL results in terms of neutral fraction (left) and temperature (right) at 1 
and 15 Myr. The third row Shows the RAMSES-RT+HLL results without the on-the-spot approximation. The fourth row shows the 
RAMSES-RT+GLF results. The bottom row shows the C^-Ray results for comparison. 
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Figure 17. 1106 test 3: RAMSES-RT+HLL versus C^-Ray comparinson. (a) Evolution of the position and speed of the I-front along the 
X-axis through the center of the box. The position plot (top) shows the x-position where xhii = 0-5, with respect to the center of the 
cloud, xc = 5 kpc, in units of the Stromgren length inside the cloud, is,cl = 0-78 kpc. The dotted horizontal lines mark the edges of the 
cloud. The speed (bottom) is plotted in units of twice the isothermal sound speed in the cloud at T = 10^ K, 2cs^i(10'^ K) = 2.35 x 10^ 
cm/s. (b) Evolution of the average ionized fraction (top) and temperature (bottom) inside the dense cloud, (c) Profiles along the x-axis 
through the box center of the ionization state (top) and temperature (bottom), at 1, 3 and 15 Myr. 
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Figure 18. 1106 test 3: RAMSES-RT+HLL versus C^-Ray comparison. Histograms of neutral fraction (top row) and temperature (bottom) 
inside the dense cloud at 1, 3 and 15 Myr (from left to right). 



tween the codes. RAMSES-RT doesn't particularly stand out, 
though, and is most similar to that of Coral at 15 Myr. The 
temperature profile for RAMSES-RT also differs notably from 
that of ATON, where the shielded region inside the cloud is 
thicker and more step-like both in the ionized fraction and 
temperature, due to the monochromatic radiation. 

Finally, Fig. |18| shows histograms of the neutral fraction 
and temperature at 1, 3 and 15 Myr for RAMSES-RT and 
C^-Ray. The comparison (also with the other codes in 1106) 



is qualitatively similar, though there is quite 
between the individual codes in these plots. 



difference 



As with the previous tests, RAMSES-RT performs well 
here and we don't really have anything out of the ordinary 
in our results. One should keep in note though that here 
we've use the non-diffusive HLL flux function, whereas in 
most cosmological simulations it would be more natural to 
use the more diffusive GLF function to have better spheri- 
cal symmetry around radiative stellar sources, which comes 
with the price of less pronounced and shorter lived shadows 



© 0000 RAS, MNRAS 000, 000-000 



24 Rosdahl et al. 




Figure 19. 1106 test 4. Maps showing slices sX z = 0.5 Lbox of the neutral fraction and temperature at times 0.05 Myr and 0.4 Myr. 
Top row shows RAMSES-RT results with physical light speed. The middle row shows the C^-Ray results (infinite light speed). The bottom 
row shows the RAMSES-RT results with one hundred times the physical light speed. 



than HLL. As discussed in AT08 though, shadows are not 
really pronounced or long-lived in a cosmological context, 
partly because radiative sources are typically spread around 
and partly because of recombination-emitted UV photons. 
Our analogue shadow experiment with the GLF intercell 
flux function gives results similar to the same experiment 
with HLL and non-OTSA: the shadow is short-lived. One 
should therefore perhaps not be overly worried that using 
moment-based RT with the GLF flux function doesn't cast 
long-lived shadows in cosmological or galactic simulations - 
these shadows are likely never long-lived anyway. 



6.5 1106 test 4: Multiple sources in a 
cosmological density field 

This test involves the propagation of ionization fronts in a 
static hydrogen-only density field taken from a cosmological 
simulation snapshot at redshift 9. The density cube is 128^ 
cells and its width is 500 /i~^ co- moving kpc (correspond- 
ing to 50 /i~^ physical kpc). The Hubble factor is /i = 0.7. 
The initial temperature is fixed at 100 K everywhere. 16 
radiative sources are picked out corresponding to the most 
massive halos in the box and these are set to radiate contin- 



uously for 0.4 Myrs, assuming a black-body spectrum with 
effective temperature T — 10^ K. The mass-dependent ra- 
diation intensity for each halo is given in a downloadable 
table (from the RT comparison project website). Our com- 
putations are run with the GLF solver and the OTSA is 
not applied. Multi-frequency is approximated with three fre- 
quency bins of Hi, Hei and Hell ionizing photons. A full light 
speed is used (i.e. /c = 1), and an analogue run is made with 
a hundred-fold light speed (/c = 100). 

Fig. 1 19| shows box slices, dX z — 0.5 Lbox, of the neutral 
fraction and temperature at times 0.05 and 0.4 Myr. Shown 
are our two runs with different light speed fractions (top 
and bottom row) , and for comparison we show the result for 
the C^-Ray code, from 110^ the I-fronts and photo-heating 
in our /c = 1 run clearly lag behind the C^-Ray result, and 
there is also less heating of the ionized gas. This is in accor- 
dance with the ATON results described in AT08, where a sim- 
ilar delay was found. They prescribed this delay to the fact 



^ Note that 1106 have likely mislabeled the maps showing the 
results from this test; their text and captions indicate the maps 
to be at 0.2 Myr, but judging from the downloadable data they 
are at 0.4 Myr. 
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Figure 20. 1106 test 4. (a) Time evolution of the mass weighted and volume weighted average ionized fractions, (b) Histograms of 
neutral fraction (top) and temperature (bottom). 



that ATON is monochromatic, but since our mult i- frequency 
approximation (three photon groups) gives results that are 
still much more similar to the ATON results than those of 
C^-Ray, especially in terms of the neutral fraction maps, we 
are inclined to blame the delay on another factor, which is 
the speed of light. Our results with the speed of light set 
to one-hundred times the physical value are shown in the 
bottom row of Fig. [19] and here the results are consider- 
ably closer to those of C^-Ray in terms of the propagation of 
heating- and I-fronts, although the maximum temperature 
in the ionized gas is still colder in comparison. All four codes 
considered in the 1106 4 test use an infinite effective speed of 
light and this may give premature fronts in the immediate 
vicinity of the sources and also further away in under-dense 
regions. Thus we are perhaps not really dealing with a delay 
in RAMSES-RT, but rather premature fronts in the 1106 codes. 
As AT08 note, we are far from reaching a static state in the 
fronts in this experiment in the run-time of 0.4 Myr and we 
should expect the different light speed runs to converge to 
similar results when static state is reached. This is further 
corroborated by our I-front light crossing time analysis from 
Section 1231 



The smaller degree of photo-heating in the ionized gas 
compared to the C^-Ray results is in line with the temper- 
ature profiles from the previous tests (e.g. Fig. 14a), and 
presumably is a consequence of the different ways multi- 
frequency is approximated. Another notable difference in 
the maps in Fig. ^] is that our fronts are smoother and less 
jagged than those in C^-Ray. This is an effect of the pho- 
ton diffusion inherent in the GLF flux function used here. 
Like AT08 we find that using HLL instead gives more jagged 
fronts. 



Fig. |20a| shows the evolution of the mass- and volume- 
weighted ionized fractions, compared for the different runs. 
The RAMSES-RT run with the physical light speed gives ion- 



ized fractions which are close (both mass- and volume- 
weighted) to the ATON ones, whereas increasing the light 
speed by a factor of hundred from the physical value gives re- 
sults closer to C^-Ray (as well as the three other codes that 
ran this test in 1106). Presumably we would converge fur- 
ther towards C^-Ray in the limit of infinite light speed, but 
computational time constraints do not allow to pursue that 
investigation. This is a further hint that the correct speed 
of light is important in the non-steady regime of ionization 
fronts. 



Finally, Fig. |20b| shows neutral fraction and tempera- 
ture histograms at three times in the test. Again there is a 
strong discrepancy between the RAMSES-RT run with /c = 1 
and C^-Ray, especially at early times, and the gap all but 
closes when fc — 100 is used instead with RAMSES-RT. There 
remains some difference though in the minimum/maximum 
temperature, being smaller /larger for C^-Ray than for our 
fc = 100 run, presumably because of our rather crude multi- 
frequency approximation. 

To summarize, there is notable discrepancy between the 
RAMSES-RT results and those presented in 1106, in that the 
RAMSES-RT ionization front lags behind, which appears to 
be due to a finite speed of light. This is corroborated to 
some degree by other papers in the literature: [Wise fc Abel| 
(|2011') use a finite light speed and seem to get results which 
are slightly lagging as well. Pawlik & Schaye ( 2008 ) specifi- 



cally do a comparison between finite and infinite light speed, 
with the finite one resulting in a delay which is substantial, 
though maybe a bit less than ours, and they do comment 
on the ionization bubbles in this test being unphysically 
large with infinite light speed methods. Other sources ap- 
pear to be in confiict with our conclusions, though: |Petkova| 
& Springel (2011a) use a finite light speed and get results 



which seem to compare well with those of C -Ray. 
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Figure 21. 1109 test 5. Maps showing slices at z = of various quantities at 100 Myrs (top panel) and 500 Myrs (lower panel). In each 
panel, the top row shows the RAMSES-RT results and the lower row shows the Capreole+C-^-Ray results for comparison. 



6.6 1109 test 5: Classical HII region expansion 

We now come to the tests described in second radiative 



transfer codes comparison paper by Iliev et al. (2009), which 
we denote as 1109. This paper provides 3 code comparison 
tests to add to those in 1106, but with the important dif- 
ference that whereas the 1106 tests are pure radiative trans- 
fer post-processing tests with fixed density fields, the tests 
in 1109 are RHD tests, i.e. with the radiative transfer di- 
rectly coupled to the gas-dynamics. Thus we now switch 
from the context of post-processing RT to hydro-coupled 
RHD. Here, the pressure buildup in photo- heated gas causes 
it to expand. Typically, the I-front is initially R-type, where 
it expands much faster than the gas response to it, which 
means RT-postprocessing is a fairly good approximation. 
The I-front then begins to slow down when it approaches 
the Stromgren radius, but gets moving again when the gas 
catches up to it, and then the front is D-type, i.e. moves 
along with the expanding gas. 

As before we compare our RAMSES-RT tests results with 
those of the grid-based short characteristics ray-tracing code 



C^-Ray ( Mellema et al.|2606 ), here coupled to the Capreole 
code, which employs a Riemann solver for the hydrodynam- 
ics. As the Capreole+C^-Ray combination performs badly on 
1109 test 6 due to instabilities, we compare also in that par- 
ticular test to C^-Ray coupled to the Eulerian TVD solver 



of Trac & Pen ( 2004 ) . The test numbers continue from the 



1106 paper, thus we now come to 1109 test 5, which concerns 
the expansion of an ionization front due to a point source 
in an initially uniform-density medium. The initial setup, 
much like that of 1106 test 2, is as follows. 



^, temperature 100 K, and ionization frac- 
^ (1109 prescribes xhi = 0). The radiative 



The box cube is Lbox = 15 kpc in width. The gas is 
hydrogen only as usual, initially homogeneous with density 
nu = 10~^ cm 
tion xhi = 10 

source is in the corner of the box and the emission rate is 
Nj = 5 X 10^^ photons s~^. The spectrum is that of a black- 
body at an effective temperature of 10^ K, and as usual we 
bin the spectrum into three groups of Hi-, Hei-, and Heii- 
ionizing photons. To keep the computational cost reasonable 
we set the light speed fraction to fc = 1/100. We use the 
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Figure 22. 1109 test 5. Radial profiles at 10, 200 and 500 Myrs, compared to the Capreole+C^-Ray results. Clockwise from top left: 
ionization fractions, pressure, temperature, Mach number, atom number density. 



GLF flux function and don't apply the OTSA, i.e photons 
are emitted from gas recombinations. The box boundaries 
are reflective at the sides touching the radiative source and 
transmissive elsewhere. The simulation time is 500 Myr. The 
base resolution of the box is 64^ cells and we apply on-the- 



fly refinement on nn and xhh gradients (see Eq. 44), so that 



the ionization front has the prescribed effective resolution of 
128^ cells. 

We first compare volume dissections at z = in the sim- 
ulation cubes at 100 and 500 Myr, for the RAMSES-RT and 

The maps show, from left 
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C -Ray results, shown in Fig. 
to right, the neutral fraction, pressure, temperature, den- 
sity and mach number, M = v/cs, where Cs = y^l.4 P/p 
is the sound speed. (Unfortunately the M output is miss- 
ing from the C^-Ray results we've downloaded.) In these 
maps, the RAMSES-RT results look very similar to those of 
C^-Ray. The xni-maps show stronger ionization immediately 
around the corner source in the C^-Ray result, and corre- 
spondingly the temperature and density maps show this 
corner gas is also hotter and more diffuse in the C^-Ray re- 
sult than in RAMSES-RT. Conversely, the photo-heating region 
is somewhat further- reaching in the RAMSES-RT result than 
in C^-Ray, as can be seen in the pressure and temperature 
maps. These small differences are likely due to the different 
approaches in approximating mult i- frequency. Notably, the 
C^-Ray maps stand out in a very similar way when com- 
pared to most of the corresponding maps from other codes 
in 1109, i.e. a stronger effect close to the radiative source but 
shorter-reaching photo-heating. 

To paint a more quantitative picture. Fig. [22] compares 
radial profiles of the same quantities (xhi, P, T, nu and 
M) for RAMSES-RT and C^-Ray at 10, 200 and 500 Myr. The 
ionization state profiles (top left) indeed show C^-Ray to ion- 
ize the gas more strongly close to the radiative source, but 
RAMSES-RT to ionize more strongly beyond the I- front. The 
I- front itself is however at very similar positions at all times. 
The pressure and temperature plots show the same thing, 
but apart from these minor differences at the extreme ends 
the shapes are very similar. The density plots show that 
C^-Ray has more has more diffuse gas close to the source as 
a result of the stronger photoheating, and also it appears 
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Figure 23. 1109 test 5: Time evolution of the ionization front, 
compared to the results from the Capreole+C'^-Ray combination. 
Upper plot shows the radius of the Stromgren sphere in units of 
5.4 kpc. The lower plot shows the speed of the front propagation. 



to have a more pronounced backflow peak around 200 Myr 
(this double peak is a temporary effect of photo-heating by 
high-energy photons beyond the I-front). The smaller back- 
flow peak in RAMSES-RT is perhaps in part a relic of on-the-fiy 
refinement, though most of the codes in 1109 actually have 
backflow peaks similarly smaller than that of C^-Ray. Un- 
fortunately we can't compare the Mach profiles directly, but 
the RAMSES-RT profiles do look very similar in shape to those 
presented in 1109 (see their Fig. 15). 

Finally, Fig. |23| shows how the position and velocity of 
the I-front (defined as where the radial average of xhii is 
equal to 0.5), for RAMSES-RT and C^-Ray. The plots for the 
two codes are virtually identical, the only noticeable differ- 
ence being a slight initial lag in the front speed. One might 
attribute this to the reduced speed of light in the RAMSES-RT 
run, but actually most other codes described in 1109 have 
a very similar lag in the initial front speed compared to 
C^-Ray. 

All in all, the RAMSES-RT results for this test com- 
pare very well with most of the codes presented in 1109. 
The RAMSES-RT result differs slightly from that of C^-Ray in 
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Figure 24. 1109 test 6: time evolution of the ionization front, 
compared to the Capreole+C^ -Ray combination. 



some aspects, most notably in the form of weaker photo- 
heating and ionization close to the radiative source and 
wider I-fronts. However, these are precisely the aspects 
where C^-Ray stands out from the other codes presented 
in 1109. 



6.7 1109 test 6: HII region expansion in a r ^ 
density profile 

This test mimics a radiative source going off in a dense cloud, 
e.g. a stellar nursery. The setup is much like that of the 
preceding test 5, the main difference being that the gas is 
here inhomogeneous, the box is much smaller, Lbox = 0.8 
kpc in width, and the radiative corner source is a hundred 
times more luminous. It radiates N^ = 5 x 10^° photons s~^, 
with a T = 10^ K blackbody spectrum (softer than test 
5), which is in RAMSES-RT approximated with three (Hi-, 
Hei-, and Hell- ionizing) photon groups, that travel at fc — 
1/100 of the speed of light. As before we use the GLF flux 
function and don't apply the OTSA. The base resolution is 
64*^ cells, but on-the-fly refinement on riH and xhh gradients 
ensures the prescribed effective resolution of 128^ cells at 
ionization and shock fronts. The initial temperature is 100 
K everywhere and the running time is 75 Myr. The dense 
cloud is centered on the corner source and is set up with a 
spherically symmetric, steeply decreasing power-law density 
profile with a small flat central core of gas number density 
no ^ 



3.2 cm and radius ro = 91.5 pc: 



nii(r) 



no if r ^ ro 

no(ro/r)^ if r ^ ro. 



(46) 



The Stromgren radius for the core density, given by 
Eq. |45[ is rs ~ 70 pc, which lies within the flat core. Thus, 
the I-front makes an initial transition from R-type to D-type 
within the core, and then may accelerate back to R-type as 
it expands into decreasingly dense gas outside the core. 

We first compare the evolution of the position and speed 
of the I-front, which is plotted in Fig. [24] for RAMSES-RT and 
the Capreole+C^ -Ray combination. The I-front moves very 
quickly (R-type) to ^ 70 pc within the first fraction of a 
Myr, stops for while and then starts to expand again with 
the flow of the gas. Both the speed and position compare well 
with C^-Ray. The initial speed in C^-Ray has an apparent 



lag which is due to under-sampling in the front positions: 
other code results which are better sampled in 1109 show 
initial speeds that are virtually identical to the RAMSES-RT 
plot. The final front position in RAMSES-RT is slightly further 
out than that of C^-Ray, though very similar to at least 
three of the codes in 1109 (Flash-HC, Licorice and RSPH). It 
also appears that the C^-Ray front is starting to accelerate 
slightly at the end, whereas the RAMSES-RT front is about to 
approach constant speed; RAMSES-RT also agrees with most 
other 1109 codes on this point. 

Fig. [25] shows the overall structure of ionization and 
the gas at 25 Myr, here with a comparison between 
RAMSES-RT (upper two rows) and TVD+C^-Ray (bottom row). 
As mentioned earlier, severe instabilities appear in the 
Capreole+C^ -Ray version of this test, so we compare to the 
more stable and symmetric test with the TVD solver. In 
addition to the default RAMSES-RT run with on-the fly AMR 
refinement, we show here in the middle row results from an 
identical RAMSES-RT run with the base resolution set to 128^ 
cells and AMR refinement turned off. There are slight spher- 
ical asymmetries appearing in the top row maps, in particu- 
lar the xhii, T and Mach maps, and the middle row maps are 
presented here to show that (the first) two of these are purely 
artifacts of on-the-fiy AMR refinement. The slightly square 
shape of the inner region in the Mach map however does not 
seem to be due to refinement and is likely rather a grid arti- 
fact which is amplified by the radially decreasing density. It 
should also be noted that the other plots produced for this 
test (I-front, Fig. 



24 



26) are abso- 



and radial profiles. Fig. 
lutely identical regardless of whether on-the-fiy refinement 
is used or the full resolution applied everywhere, suggesting 
that AMR refinement produces very robust results. 

As usual the I-front is considerably wider in RAMSES-RT 
than in the C^-Ray results, though we don't find the same 
discrepancy as in the previous test between the photoheat- 
ing intensity close to the source (also, there is no such dis- 
crepancy here between C^-Ray and the other codes in 1109). 
The two maps furthest to the right, of density and Mach 
number, show the expanding shell of dense gas due to pho- 
toheating. Here the shell appears considerably thinner in 
RAMSES-RT than in TVD+C^-Ray, and indeed TVD+C^-Ray ap- 
pears to have the thickest density shell of any of the codes in 
1109 (Capreole+C^-Ray included, but here there are also se- 
vere asymmetries). The RAMSES-RT maps compare well with 
the C^-Ray ones, and to most of the maps in 1109, and don't 
show any I-front instabilities that seem to have a tendency 
to come up in this test (and 1109 do show that these are 
numerical and not physical instabilities). 

Fig. [26] shows a comparison between RAMSES-RT and 
TVD+C^-Ray for radially averaged profiles at 3, 10 and 25 
Myr of the ionization state, pressure, temperature, density 
and Mach number. The comparison is generally very good. 
The I-front (and corresponding density shock) lag a little 
behind in C^-Ray, but it actually lags a little behind all but 
one code in this test in 1109, and RAMSES-RT is spot-on com- 
pared with those others in every respect. 

All in afl, RAMSES-RT thus performs well on this test, 
and no problems appear that are worth mentioning. 
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Figure 25. 1109 test 6. Maps showing slices at z = of various quantities at 25 Myrs. The top row shows the RAMSES-RT results with 
adaptive refinement. The middle row shows results also from RAMSES-RT, but with a fully refined box and adaptive refinement turned off. 
The bottom row shows the TVD+C^-Ray results for comparison. 
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Figure 26. 1109 test 6. Radial profiles at 3, 10 and 25 Myrs, compared to the TVD+C^-Ray results. Clockwise from top left: ionization 
fractions, pressure, temperature, Mach number, atom number density. 
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6.8 1109 test 7: Photo-evaporation of a dense 
clump 

The setup of this test is identical to test 3 in 1106, where UV 
radiation is cast on a gas cloud, creating a shadow behind 
it and a slowly-moving I- front inside it. Here however, since 
the hydrodynamics are turned on, photo-heating causes the 
cloud to expand outwards and simultaneously contract at 
the center. We recap the setup: 

The box is Lbox = 6.6 kpc in width. A spherical cloud 
of gas with radius rdoud = 0.8 kpc is placed at (xc, 2/c, Zc) — 
(5,3.3,3.3) kpc from the box corner. The density and tem- 
perature are n^^ = 2 x 10"^ cm"^ and T°^* = 8000 K 
outside the cloud and ng°^^ = 200 ng^' = 4 x 10"^ cm-^ 
and T'^^^^'^ = 40 K inside it. From the x = boundary 
a constant ionizing flux oi F — 10^ photons s~^ cm~^ 
is emitted towards the cloud. The spectrum is a black- 
body with Teff = 10^ K, which is approximated with three 
(Hi— , Hel—, Hell— ionizing) photon groups, that travel at 
fc — 1/100 of the speed of light. The simulation time is 
50 Myr, considerably longer than the 15 Myr in the corre- 
sponding pure RT test. The base resolution is 64^ cells, but 
on-the-fly refinement on nn, x^i and xhh gradients ensures 
the prescribed effective resolution of 128^ cells at ionization 
and shock fronts. In order to best capture the formation of 
a shadow behind the cloud, we focus on a RAMSES-RT run 
with the HLL solver, but we also show some results with 
the usual GLF solver. The OTSA is applied in this test. 

Fig. [27| shows slices in the x^-plane through the mid- 
dle of the box of various quantities at 10 and 50 Myr, for 
the RAMSES-RT result and C^-Ray for comparison. As in the 
corresponding pure RT test, it can be seen from the xhi 
maps that the shadow behind the cloud less conserved with 
RAMSES-RT than with C^-Ray, though the HLL solver does 
a much better job though than GLF. However, the dif- 
fusion of photons doesn't have a large impact on the re- 
sulting dynamics, or even the propagation of the Lfront 
along the axis of symmetry. The shadow becomes thinner 
towards the end of the run with all codes in 1109, though 
it is thinner than most in RAMSES-RT+HLL, and it pretty 
much disappears in RAMSES-RT+GLF. The shadow thickness 
in RAMSES-RT+HLL is still comparable at 50 Myrs to the re- 
sults of RSPH, Zeus-MP and Licorice in 1109. The pressure 
maps of RAMSES-RT+HLL, C^-Ray and other codes in 1109 are 
very similar both at 10 and 50 Myrs, though C^-Ray, and 
also to some extent Flash-HC and Licorice have a fork- 
like shape inside what remains of the shadow at 50 Myr. 
The other codes have the same shape as RAMSES-RT+HLL 
in this region. The temperature maps are similar as well, 
though there appears here to be slightly more heating in 
RAMSES-RT both inside the cloud and in the expanding den- 
sity shell around it. This is because for some reason the 
initial density, both inside and outside the cloud, is slightly 
too high in the C^-Ray run, which also seems to result in 
slightly too low initial temperatures. The shell expands in 
a very similar way for the two codes, as can be seen in the 
density and Mach slices. The expansion goes a bit further, 
though, in RAMSES-RT. Also, the expanding cloud seems to 
develop a slightly hexagonal shape in RAMSES-RT, an effect 
which is not apparent in any of the codes in this test in 
1109 (though there is a hint of it in the Flash-HC result). 
It can only be speculated that this is a grid artifact. To 



be sure it doesn't have to do with the on-the-fly refine- 
ment we ran an identical experiment with a base resolution 
of 128^ cells and no refinement in RAMSES-RT+HLL. The 
RAMSES-RT+HLL maps and plots presented here are virtu- 
ally identical to this non-refinement run, except of course for 
graininess in the slice maps. None of these discussed effects 
(hexagons and a slightly over-extended I-front compared to 
other codes) are thus due to on-the-fly refinement. 

Next we turn our attention to the evolution of the po- 
sition and speed of the I-front along the x-axis of symmetry 
through the box. This is presented f or th e RAMSES-RT (HLL 
and GLF) and C^-Ray runs in Fig. 28a The I-front prop- 



agation is considerably different between RAMSES-RT and 
C^-Ray, but actually C^-Ray considerably stands out from 
other codes in this test in 1109. For the first 7 Myrs or so, 
the RAMSES-RT front lags behind that of C^-Ray and in fact 
all the codes in 1109. This is due to the reduced speed of light: 
before hitting the cloud, the photons have to travel from the 
left edge of the box through a very diffuse medium - so dif- 
fuse that here the I-front speed apparently is approaching 
the speed of light, or is at least considerably faster than the 
one-one-hundredth of the light speed which is used in the 
RAMSES-RT run. However, once the I-front in the RAMSES-RT 
run has caught up, the reduced light speed should have a 
negligible effect on the results. After roughly 7 Myr, the 
RAMSES-RT I-front overtakes C^-Ray front, and stays ahead 
of it for the remainder of the run. This however is also the 
case for most of the codes in 1109; their I-front is ahead of 
the C^-Ray front, and four out of six codes end up with the 
Lfront at - 5.6 kpc. The RAMSES-RT+HLL front ends up 
at ^ 5.7 kpc, so slightly ahead of what is typically found 
in 1109. Using the GLF solver instead of HLL has the effect 
that the I-front disappears soon after 40 Myrs, which is due 
to diffusive photons eating into the shadow from it's edges, 
but up to that point the I-front evolution is much the same. 
RAMSES-RT also reproduces the retreat of the I-front between 
roughly 30 and 40 Myrs, which is seen in all runs in 1109. 
This momentary negative speed is due to the expansion of 
the cloud and the D-type movement of the I-front with the 
gas. 

Fig. |28b| shows histograms of the gas temperature and 
Mach number at 10 and 50 Myr in the RAMSES-RT+HLL 
and C^-Ray runs. The shapes of the histograms are very 
similar between the two codes (and are also very similar 
to RAMSES-RT+GLF, which is not shown). The temperature 
histograms are a bit offset, but this is also the case when 
the C^-Ray histogram is compared with any of the others 
in 1109, and is probably just due to the slightly off starting 
temperature in the C^-Ray initial conditions. 

Finally, Fig. [29l shows a comparison between RAMSES-RT 
and C^-Ray of profiles along the x-axis of symmetry of the 
various quantities at 1, 10 and 50 Myrs. The profiles compare 
badly at 1 Myr, but as already discussed this is simply due 
to the I-front having not caught up at this early time when 
using the reduced speed of light. At later times the profiles 
generally compare well, though we see these effects which 
have already been discussed, of a further expanding density- 
front out of the original cloud, a further progressed I-front, 
and slight differences in density and temperature which are 
due to slightly wrong initial conditions in the C^-Ray run. 
The RAMSES-RT profile plots show a staircase effect which 
is most obvious in the 50 Myrs plot at the radial interval 
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Figure 27. 1109 test 7. Maps showing slices at z = 0.5 Lbox of various quantities at 10 Myrs (top panel) and 50 Myrs (lower panel). 
In each panel, the top row shows the RAMSES-RT+HLL results, the middle row shows RAMSES-RT+GLF and the bottom row shows the 
Capreole+C^-Ray results. 
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Figure 28. 1109 test 7. (a) Time evolution of the position (top) and speed (bottom) of the ionization front along the x-axis of symmetry 
through the center of the box. (b) Histograms of the gas temperature (upper panel) and flow Mach number (lower panel) at 10 and 50 
Myr for RAMSES-RT and Capreole+C^-Ray. 
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Figure 29. 1109 test 7. Proflles along the x-axis of symmetry through the center of the box, at 1, 10 and 50 Myr for the RAMSES-RT and 
Capreole+C^-Ray results. Clockwise from top left: ionization fractions, pressure, temperature, Mach number, atom number density. 



0.45 < r/Lbox ^ 0.75: this is simply due to the grid being 
unrefined at this x-interval along the axis of symmetry, i.e. 
at the effective base resolution of 64^ cells per box width. 
The run with the full resolution and no AMR refinement 
shows no staircases, but otherwise the results are identical 
to those shown here. 



We have made an alternative run with RAMSES-RT+HLL 
with the speed of light fraction set to fc = 1/10 rather than 
the default 1/100, and here the initial evolution of the I-front 
position and radial profiles at 1 Myr are almost identical to 
those of C^-Ray. At later times the results are very much in 
line with those where fc — 1/100, except the I-front position 



is slightly more advanced at 50 Myr, or at 5.78 kpc rather 
than at 5.71 kpc. 



In summary, RAMSES-RT performs well on this test with 
no apparent problems. The reduced light speed (/c = 1/100) 
has very little effect on the results and on-the-fly refinement 
gives results which are identical to the fully refined simu- 
lation with a homogeneous 128^ cells grid. Even using the 
diffusive GLF solver retains much of the results (I-front de- 
velopment, cloud expansion), except that the I-front disap- 
pears a bit prematurely. 
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6.9 RT test conclusions 

RAMSES-RT performs very well on all the tests from 1106 and 
1109, with no discrepancies to speak of from expected results 
or those from other codes. The most notable discrepancies 
clearly result from the reduced speed of light approximation 
used in RAMSES-RT, though it is not always clear whether this 
gives results that are 'less correct' than when using infinite 
speed of light, as is the case for most of the codes compared 
in the RT comparison project. Our shadows are also consid- 
erably shorter lived with the GLF intercell flux function than 
those of the RT comparison project codes (most of which use 
ray-tracing schemes) . This can be fixed for problems involv- 
ing shadows and idealized geometries by using the HLL flux 
function instead, but as we showed in §3.2| the sacrifice is 
that isotropic sources become not so isotropic. Many codes 
in the RT comparison project show various instabilities and 
asymmetries in ionization fronts; no such features are man- 
ifested in the RAMSES-RT results. 



7 DISCUSSION 

In this paper, we have presented a new implementation of 
radiation hydrodynamics in the RAMSES code. It is based on a 
moment representation of the radiation field, where we have 
used the Ml closure relation to define a purely local vari- 
able Eddington tensor. Because the resulting system is a set 
of hyperbolic conservation laws, we have exploited the Go- 
dunov methodology to design a time-explicit, strictly photon 
conserving radiation transport scheme. The resulting algo- 
rithm is first order accurate in space and time, and uses var- 
ious Riemann solvers (GLF and HLL) to compute radiation 
fluxes. The main novelty compared to our previous imple- 
mentation (AT08) is the coupling between gas and radia- 
tion, resulting in a fully consistent radiation hydrodynamics 
solver, and the introduction of adaptive mesh techniques in 
the radiation transport step, making use of both the AMR 
and parallel computing capabilities of RAMSES. Overall, the 
code was quite easy to implement, owing to the explicit na- 
ture of the time integration scheme. The price to pay is the 
need to resolve the propagation of hyperbolic waves travel- 
ing at or close to the speed of light. Among many different 
options available to overcome this constraint, we have chosen 
to use the "reduced speed of light" approximation. This ap- 
proximation is valid when the propagation speed of I-fronts 
is still slower than the (reduced) light speed. We have devel- 
oped a recipe to assess the validity of this approximation, 
based on the light-crossing time of Stromgren spheres. We 
have verified that this framework indeed allows us to esti- 
mate in advance the speed-of-light reduction factor reliably. 
We have shown, for example, that in cosmological problems, 
such as cosmic reionization, using the correct value for the 
speed of light is crucial, and using either a reduced or an 
infinite speed of light (like in some ray-tracing codes) might 
result in large inaccuracies. 

This new algorithm has already been used in galaxy 
formation studies, exploiting the coupling between radia- 
tion and hydrodynamics offered by RAMSES-RT. In Rosdahl 



& Blaizot ( 2012 ), we studied the impact of ionizing radiation 



in determining the thermal state of cold filaments stream- 
ing into high redshift galaxies, allowing us to make accu- 
rate observational predictions and demonstrating a possible 



link between cold streams and Lyman- alpha blobs. More re- 
cently, we have also explored the role of ionizing radiation 
in the overall efficiency of stellar feedback (Geen et al. 2013, 
Powell et al. 2013, both in preparation). Beyond ionizing ra- 
diation, possible extensions of RAMSES-RT are the inclusion 
of photo-dissociating radiation and the thermochemistry of 
molecules, as well as the effect of dust as an additional source 
of opacity and thermal regulation inside star forming galax- 
ies. This would require introducing additional photon groups 
(such as a far UV and IR photons) and the associated mi- 
crophysics, but the overall methodology would remain very 
similar. 

In order to improve the current algorithm, we have 
many possibilities ahead of us. One obvious development 
is to develop a second-order sequel of our current first order 
Godunov solver. Second-order Godunov schemes, both in 
time and space, are used routinely in hydrodynamics codes 
(such as the MUSCL scheme in RAMSES). This might reduce 
significantly the rather large diffusivity of our current im- 
plementation. However, since photo-ionization and photo- 
dissociation problems are governed to a large extent by the 
thermo-chemistry, it is not clear how much the accuracy of 
the results would depend on the advection scheme. A second 
route we would like to explore in the future is the optional 
introduction of radiation sub-cycles during each adaptive hy- 
dro step. This is quite challenging since it would in principle 
require decoupling in time of the various AMR levels, result- 
ing in the loss of strict photon conservation. In some cases, 
however, it is advantageous to sacrifice the exact number 
conservation of photons in favour of modelling the correct 
speed of light with many radiation sub-cycles. In any case, 
this would offer us a new tool with greater fiexibility. Along 
the same lines, because of the fundamentally different prop- 
agation properties of I-fronts in the IGM on one hand and 
deep inside galaxies on the other, we could couple RAMSES-RT 
to ATON: use ATON to transport radiation on the coarse grid 
with GPUs at the full speed of light, and use RAMSES-RT on 
the fine AMR levels at a reduced light speed. This would 
require us to define two photon group populations that mir- 
ror each other: a large scale, low density photon population 
that propagates at the correct speed of light and makes use 
of GPU acceleration (if available), and a small scale, high 
density photon population that makes use of the "reduced 
speed of light" approximation. Coupling properly the two 
photon group populations will of course be quite challeng- 
ing and at the heart of this new avenue of research. A last 
development we have in mind is the introduction of radia- 
tion pressure as a new channel of coupling radiation with 
hydrodynamics. This is highly relevant for studies focusing 
on radiation pressure on dust, from both young star clusters 
and supermassive black holes. 



ACKNOWLEDGEMENTS 

We are grateful for the help and insight provided by 
Stephanie Courty, Julien Devriendt, Yohan Dubois, and 
Leo Michel-Dansac. This work was funded in part by the 
Marie Curie Initial Training Network ELIXIR of the Eu- 
ropean Commission under contract PITN-GA-2008-214227, 
by the European Research Council under the European 
Unions Seventh Framework Programme (FP7/2007-2013) 



© 0000 RAS, MNRAS 000, 000-000 



34 Rosdahl et al. 



I ERC Grant agreement 278594-GasAroundGalaxies, and 
the Marie Curie Training Network CosmoComp (PITN- 
GA-2009-238356). The tests were performed using the 
HPC resources of CINES under the allocation 2011- 
c201 1046642 made by GENCI (Grand Equipement National 
de Calcul Intensif), and the Cray XT-5 cluster at CSCS, 
Manno, Switzerland. We also acknowledge computing re- 
sources at the CC-IN2P3 Computing Center (Lyon/Villeur- 
banne - France), a partnership between CNRS/IN2P3 and 
CEA/DSM/Irfu. JB acknowledges support from the ANR 
BINGO project (ANR-08-BLAN-0316-01). 



REFERENCES 

Abel T., Norman M. L., Madau P., 1999, ApJ, 523, 66 
Abel T., Wandelt B. D., 2002, MNRAS, 330, L53 
Agertz O. et al., 2007, Monthly Notices of the Royal As- 
tronomical Society, 380, 963 
Altay G., Croft R. A. C, Pelupessy I., 2008, MNRAS, 386, 

1931 
Alvarez M. A., Bromm V., Shapiro P. R., 2006, ApJ, 639, 

621 
Anninos P., Zhang Y., Abel T., Norman M. L., 1997, New 

Astronomy, 2, 209 
Aubert D., Teyssier R., 2008, MNRAS, 387, 295 
Aubert D., Teyssier R., 2010, ApJ, 724, 244 
Back S., Di Matteo P., Semelin B., Combes P., Revaz Y., 

2009, A&A, 495, 389 

Back S., Semelin B., Di Matteo P., Revaz Y., Combes P., 

2010, A&A, 523, 4 

Black J. H., 1981, MNRAS, 197, 553 

Bruzual G., Chariot S., 2003, MNRAS, 344, 1000 

Cantalupo S., Porciani C, 2011, Monthly Notices of the 

Royal Astronomical Society, 411, 1678 
Cantalupo S., Porciani C, Lilly S. J., Miniati P., 2005, 

ApJ, 628, 61 
Gen R., 1992, ApJS, 78, 341 
Gen R., 2002, AJS, 141, 211 
Gen R., Fang T., 2006, ApJ, 650, 573 

Giardi B., Ferrara A., Marri S., Raimondo G., 2001, MN- 
RAS, 324, 381 
Giardi B., Stoehr P., White S. D. M., 2003, Monthly Notice 

of the Royal Astronomical Society, 343, 1101 
Commercon B., Teyssier R., Audit E., Hennebelle P., 

Chabrier G., 2011, A&A, 529, 35 
Croft R. A. C, Altay G., 2008, MNRAS, 388, 1501 
Dijkstra M., Loeb A., 2009, MNRAS, 396, 377 
Dubroca B., Feugeas J., 1999, Comptes Rendus de 

r Academic des Sciences - Series I - Mathematics, 329, 

915 
Finlator K., Ozel P., Dave R., 2009, MNRAS, 393, 1090 
Gnedin N. Y., 2000, ApJ, 535, 530 
Gnedin N. Y., Abel T., 2001, New Astronomy, 6, 437 
Gnedin N. Y., Ostriker J. P., 1997, Astrophysical Journal 

V.486, 486, 581 
Gonzalez M., Audit E., Huynh P., 2007, A&A, 464, 429 
Haiman Z., Thoul A. A., Loeb A., 1996, Astrophysical 

Journal v.464, 464, 523 
Harten A., Lax P. D., Leer B. v., 1983, SIAM review, 25, 

35 



Hopkins P. P., Keres D., Murray N., Quataert E., Hern- 

quist L., 2012, Monthly Notices of the Royal Astronomical 

Society, 427, 968 
Hopkins P. P., Quataert E., Murray N., 2011, MNRAS, 

417, 950 
Hui L., Gnedin N. Y., 1997, MNRAS, 292, 27 
Iliev L T. et al., 2006a, MNRAS, 371, 1057 
Iliev L T., Mellema G., Pen U.-L., Merz H., Shapiro P. R., 

Alvarez M. A., 2006b, MNRAS, 369, 1625 
Iliev I. T. et al., 2009, MNRAS, 400, 1283 
Krumholz M. R., Klein R. I., McKee C. P., 2012, ApJ, 754, 

71 
Krumholz M. R., Klein R. I., McKee C. P., Bolstad J., 

2007, ApJ, 667, 626 
Laursen P., Sommer-Larsen J., 2007, ApJ, 657, L69 
Leitherer C. et al., 1999, AJS, 123, 3 

Levermore C. D., 1984, Journal of Quantitative Spec- 
troscopy and Radiative Transfer, 31, 149 
Martel H., Shapiro P. R., 1998, MNRAS, 297, 467 
MaselU A., Giardi B., Kanekar A., 2009, MNRAS, 393, 171 
MaselU A., Ferrara A., Giardi B., 2003, MNRAS, 345, 379 
Mellema G., Iliev I. T., Alvarez M. A., Shapiro P. R., 2006, 

New Astronomy, 11, 374 
Mihalas D., Mihalas B. W., 1984, New York 
Miralda-Escude J., Haehnelt M., Rees M. J., 2000, ApJ, 

530, 1 
Nakamoto T., Umemura M., Susa H., 2001, MNRAS, 321, 

593 
Ocvirk P., Aubert D., 2011, MNRAS, 417, L93 
Oppenheimer B. D., Schaye J., 2013, arXiv, 1303, 19 
Osterbrock D. E., Ferland G. J., 2006, Astrophysics of 

gaseous nebulae and active galactic nuclei 
Pawlik A. H., Schaye J., 2008, MNRAS, 389, 651 
Pawlik A. H., Schaye J., 2009, MNRAS, 396, L46 
Pawlik A. H., Schaye J., 2011, Monthly Notices of the Royal 

Astronomical Society, 412, 1943 
Petkova M., Springel V., 2009, Monthly Notices of the 

Royal Astronomical Society, 396, 1383 
Petkova M., Springel V., 2011a, Monthly Notices of the 

Royal Astronomical Society, 415, 3731 
Petkova M., Springel V., 2011b, Monthly Notices of the 

Royal Astronomical Society, 412, 935 
Pierleoni M., Maselh A., Giardi B., 2009, MNRAS, 393, 

872 
Powell L. C., Slyz A., Devriendt J., 2011, MNRAS, 414, 

3671 
Razoumov A. O., Cardall C. Y., 2005, MNRAS, 362, 1413 
Reynolds D. R., Hayes J. C, Paschos P., Norman M. L., 

2009, Journal of Computational Physics, 228, 6833 
Rijkhorst E.-J., Plewa T., Dubey A., Mellema G., 2006, 

A&A, 452, 907 
Rosdahl J., Blaizot J., 2012, MNRAS, 2837 
Scannapieco C. et al., 2012, Monthly Notices of the Royal 

Astronomical Society, 423, 1726 
Shapiro P. R., Iliev I. T., Alvarez M. A., Scannapieco E., 

2006, ApJ, 648, 922 
Sokasian A., Yoshida N., Abel T., Hernquist L., Springel 

v., 2004, MNRAS, 350, 47 
Sun P., 1996, Journal of Computational Physics, 127, 152 
Susa H., 2006, Publications of the Astronomical Society of 

Japan, 58, 445 
Tasker E. J., Brunino R., Mitchell N. L., Michielsen D., 



© 0000 RAS, MNRAS 000, 000-000 



RAMSES-RT 35 



Hopton S., Pearce F. R., Bryan G. L., Theuns T., 2008, 

Monthly Notices of the Royal Astronomical Society, 390, 

1267 
Teyssier R., 2002, A&A, 385, 337 
Toro E. F., 1999, Riemann Solvers and Numerical Methods 

for Fluid Dynamics: A Practical Introduction. Berlin 
Trac H., Pen U.-L., 2004, New Astronomy, 9, 443 
Vaytet N., Audit E., Dubroca B., 2010, Numerical Model- 
ing of Space Plasma Flows, 429, 160 
Verhamme A., Schaerer D., Maselh A., 2006, A&A, 460, 

397 
Verner D. A., Ferland G. J., Korista K. T., Yakovlev D. G., 

1996, Astrophysical Journal v.465, 465, 487 
Wadsley J. W., VeeravalU G., Couchman H. M. P., 2008, 

Monthly Notices of the Royal Astronomical Society, 387, 

427 
Whalen D., Norman M. L., 2006, AJS, 162, 281 
Wise J. H., Abel T., 2008, ApJ, 685, 40 
Wise J. H., Abel T., 2011, MNRAS, 414, 3458 
Yajima H., Li Y., Zhu Q., Abel T., 2012, Monthly Notices 

of the Royal Astronomical Society, 424, 884 
Zahn O., Lidz A., McQuinn M., Dutta S., Hernquist L., 

Zaldarriaga M., Furlanetto S. R., 2007, ApJ, 654, 12 



APPENDIX A: RAMSESRT 
NON-EQUILIBRIUM THERMOCHEMISTRY 

We describe here in detail the non-equilibrium thermo- 
chemistry we have implemented for RAMSES-RT to accom- 
modate for the interactions between photons and gas. 
A thermochemistry step in RAMSES-RT considers a sin- 
gle cell of gas at a time with a given state hi — 
{p,pu,E,pxHiupxHeiupxHeui,Ni,Fi) (respectively, mass 
density, momentum density, energy density, hydrogen and 
helium ion abundances, photon densities and fluxeaj) and 
evolves numerically over a time-step At the thermochemical 
state 

VIt = (^, a^Hii, a^Hell, a^Heiii, A^i , •• A^M , i^i, •• Fm) 

(Al) 
(where e — E — l/2p\j? is the thermal energy density), i.e. 
solves the set of 4 + 2M coupled equations 



dUr 
dt 



S, 



(A2) 



where S = Ut ■ 

Due to the stiffness of the thermochemistry equations, 
it is feasible to solve them implicitly, i.e. using S{U^^^) on 
the right hand side (RHS), which guarantees stability and 
convergence of the solver. However a fully implicit solver is 
complicated in implementation, computationally expensive 
and not easily adaptable to changes, e.g. a varying number 
of photon groups or additional ion/chemical abundances. In- 
stead we take an approach inspired by |Anninos et al. ( 1997). 
The idea is to solve one equation at a time in a specific order. 



" Here we ignore the metal mass density, which is optionally 
stored in every cell, but at this time is not used in the non- 
equilibrium thermochemistry. 



and on the RHS use forward-in-time (FW) values, i.e. evalu- 
ated at t + At, wherever available, but otherwise backwards- 
in-time (BW) values, evaluated at t. So for the first variable 
we choose to advance in time, there are no FW variables 
available. For the next one, we can use the FW state of the 
first variable, and so on. In that sense the method can be 
thought of as being partially implicit. 

The cell-thermochemistry is called once every RT time- 
step of length Ati^T, but in each cell it is split into local 
sub-steps of length At that adhere to the 10% rule. 



AUt 



Ut 



^0.1, 



(A3) 



where AUt is the change in Ut during the sub- step. The 
RT step thus contains a loop for each cell, which calls the 
thermo_step(^T, At) routine once or more often: first with 
At = AtRT, then possibly again a number of times to fill 
in AtRT if the first guess at At proves too long to meet the 



condition set by ( A3 ) 



The thermo_step(Z//T, At) routine performs the following 



tasks: 



i) N and F update 

ii) E update 

iii) xhi update 

iv) XHeii and XHeiii update 

v) Check if we are safe to use a bigger time-step 



Tasks (ii) to (iv) are in the same order as in Anninos 



et al. (1997), but they don't include radiative transfer in 



their code, so there is no photon update. The argument we 
have for putting it first rather than anywhere else is that the 
photon densities appear to be the most dynamic variables 
and so are also most likely to break the time-step condition 



( A3 ) . This we want to catch early on in the thermochemistry 
step so we avoid doing calculations of tasks (ii) to (v) that 
turn out to be useless because of the too-long time-step. 

We now describe the individual tasks. Temperature de- 
pendent interaction rates frequently appear in the tasks - 
their expressions are given in Appendix^ The temperature 
can at any point be extracted from the energy density and 
ionization state of the gas via 



pkB 



(A4) 



where 7 is the ratio of specific heats (usually given the value 
of 5/3 in RAMSES, corresponding to monatomic gas), uih 
the proton mass, fe the Boltzmann constant and p is the 
average mass per particle in the gas, in units of tuh- 



|(i)| Photon density and flux update 

The photon number densities and fluxes, Ni and Fi, are up- 
dated one photon group i at a time. For the photon density 
the equations to solve are 



dt 



= N^+C^-N^ Di 



(A5) 



where Ni represents the time derivative of Ni given by 
the RT transport solver (which is nonzero only if the 
smoothed RT option is used), d represents photon-creating 
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re-combinations, and Di represents photon-destroying ab- 
sorptions. The creation term is non-existent if the OTSA is 
used (emitted photons are assumed to be immediately reab- 
sorbed), but is otherwise given by 



HlI.Hell.Helll 



a 



E 



/ A Bn 

(ttj — aj ) rij Tie 



(A6) 



where the 6J|^ factor is a boolean (1/0) that states which 
photon group j-species recombinations emit into and a^ 
and Q^ are the temperature dependent case A and B re- 
combination rates for the recombining species. The photon 
destruction factor is given by 



Hl,HeI,HeII 



D,, 



E 



(A7) 



where Cr is the (reduced) light speed and afj is the cross- 
section between species j and photons in group i. 

Photon emission from recombination is assumed to be 
spherically symmetric, i.e. to go in all directions. It is there- 
fore purely a diffusive term, and the photon flux equation 
only includes the photo- absorbtions: 

dt 



F,D^ 



(A8) 



where Fi is the time derivative used only in smoothed RT 



and the destruction factor remains as in (A7) 



Equations ( A5 ) and ( A8 ) are solved numerically using 



a partly semi-implicit Euler (SIE) formulation, in the sense 
that they are semi-implicit in the photon density and flux 
but otherwise explicit (in temperature and the ion abun- 
dances). A tiny bit of algebra gives: 



N^ 



t-\-At 



:it-\-At 



AT/ + AtjN, + a] 
1 + AtDi 

F/ + Am 



(A9) 



(AlO) 



1 + AtA 

where all the variables at the RHS are evaluated at the be- 
ginning of the time-step, i.e. at t. 

For each photon group update, the 10% rule is checked: 
if 

' ' ^. ^' >0-l. (All) 

the cool_step routine returns with an un-updated state 
but instead a recommendation for a new time-step length 
At new = 0.5 At, SO the routine can be called again with a 
better chance of completing. 



(A12) 



(ii) Thermal update 

Due to the dependency of /i on the ionization fractions it is 
easiest to evolve the quantity 

T 
T = — 

where fi can be extracted via 

M = [X(l + XHll) + 1^/4(1 + XHell + 2XHeIIl) ]"' , (A13) 

with X and Y — 1 — X the hydrogen and helium mass 
fractions, respectively. Here we ignore the metal contribution 
to /i, which in most astrophysical contexts is negligible. 



The temperature is updated by solving 



dTa 



(7 - l)mH 



A, 



(A14) 



dt pkB 

where A = e = 7i-\-C,7iis the photoheating rate and C the 
cooling rate. These rates are calculated as follows: 

The photoheating rate Ti is a. sum of the heating contri- 
butions from all photoionization events: 



Hl,Hel,Hell 

3 



Hel,Hell „qq 

V nj aj{v)F{v)[hv-t,]dv, (A15) 



where v is photon frequency, F{h') local photon flux and 
Cj photoionization energies. With the discretization into M 



(AI6) 



photon groups, (A15) becomes 

Hl,HeI,HeII M 

3 *=1 

where e^, afj^ afj and are the photon average energies, av- 
erage cross sections and energy weighted cross sections, re- 
spectively, for ionization events between group i and species 
j (see Eqs. [9pT ). 

The primordial cooling rate C is given by 

C = [Chi(T) + ^Hi(T)] Tie nm (A17) 

+ CHel(T) Tie TlHel 

+ CHell(T) + '0HeIl(T) + ?7HeII (^) + ^Hell (T) neUUeU 

+ mii{T) rie nun 

+ mem(T) Tie nueiu 

+ 0{T) He (nun + nHell + 4nHeIIl) 

+ UJ{T) Tie, 

where the various cooling processes are collisional ioniza- 
tions ^5 collisional excitations ip, recombinations 77, dielec- 
tronic recombinations cj, bremsstrahlung and Compton 
cooling I^7, all analytic (fitted) functions of temperature 
taken from various sources. The complete expressions are 
listed (with references) in Appendix [E| If the OTSA is used, 
the T]^ coefficients are replaced with 77^. 



The temperature update (A14) is solved numerically 



using semi-implicit formulation in T^, using FW values of 
photon densities and BW values of H and He species abun- 
dances. The temperature is updated to 

^^^* (AI8) 



^ ~^^l-A^i^At' 

where K = ^ • The temperature-derivative, A^ = ^-, 
is found by algebraically differentiating each of the primor- 
dial cooling rate expressions in the case of C (and using 
■§^ — IJL^). The temperature derivative of the heating rate 
is zero. 

With T^^^* in hand, the time-stepping condition is 
checked, i.e if 



|y^+At 



rpt I 



rrt 



>0.1, 



(A19) 



cool_step is re-started with half the time-step length. In 
tests we've found that the usual time-step constraint given 
here is not enough to ensure stability, as the temperature in 
some cases oscillates, even in a divergent way. A and A^ are 
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both evaluated backwards in time, i.e. at t, and the large 
difference that can exist in these values from t to t + At 
appears to cause these instabilities. To fix that we include 
also a first-order time- step constraint on the temperature, 
i.e. if 



|KAAt| 



>0.1, 



(A20) 



the time-step length is halved. With this fix, we have not 
seen further temperature oscillations, but there is no guar- 
antee that numerical instabilities are eliminated. 



I(iii) I Hydrogen ionized fraction update 

The Hii abundance is affected by collisional ionizations, pho- 
toionizations, and recombinations, i.e. 



dnmi 
dt 



riHi I /^mrie + '^cFmiCrNi 1 - nHiittmi^e, (A21) 



where /3hi (T) is the rate of collisional ionizations by electrons 
and (y-mi{T) the case A hydrogen recombination rate, which 
is replaced here by aHii if ^h^ OTSA is used. In terms of 



ionization fraction, (A21) becomes 

M 



dXY 



dt 



(1 - XHll) 



Pmrie + ^ afu^CrNi 



(1 - XHll) C - XHII D 

C-xuu {C + D), 



XHllttHlI^e 



(A22) 



where we have in the second line separated the rates into Hii 
creation C and destruction D, and in the third line collected 
multiples of xhii- 

To prevent stiffness-induced instabilities, we have gone 
for an approach which is semi-implicit in xhii^ 



where 



J - 



-^HII 



dxnii 

dxmi 

dC 



: Xhii + At 



C-x\,,,{C + D) 



{C + D)- xuu 



J At 



dC dD 

+ 



(A23) 



(A24) 



and the creation and destruction derivatives are given by 



dC 
dxnu 

dP 

dxYLll 






^HttHii - rieT^ii X 



2^^ttHlI 



dT ■ 



(A25) 
(A26) 



We end with the usual check if the 10% rule is broken. 



i.e. if 



^t+At 



a^Hi: 



>0.1, 



(A27) 



cool_step is restarted with half the time-step length. Like 
with the temperature a first-order check is also made, i.e. 



\C - x\,,,{C ^ D)\ 



At > 0.1. 



(A28) 



(iv) Helium ionized fractions update 



Though the Hei fraction is not a cell variable (it can be 
obtained via XHei = 1— a^Heii— a^Heiii), it is evolved in order to 
make a consistency check at the end of the helium updates. 
Before each of the helium fraction updates, we recalculate 
Ue and /i to reflect the new FW abundances. 
The Hei fraction is set by 



dx^ 



dt 



A 
a^HellQ^He 



XUel I /^Hel^e " ^ af^^^CrNi 



C - XHel D, 



(A29) 



i.e. Hell recombinations and collisional- and photo- 
ionizations of Hei. As usual a^ is replaced by a^ in the case 



of the OTSA. In the second line of (A29) we've separated 



the RHS into Hei creation C and destruction D. 

Here we follow Anninos et al. (1997*) and do the Hei 
update with 



t+At _ ^Hel + CAt 



(A30) 



The update is partly implicit, since it uses updated values 
of Ar*+^*, T^+^*, and x^+^* (^ /i and rie), but un-updated 
values of Xneii and XHem- 

We then evolve the Hell fraction. The differential equa- 
tion to solve is 



dxi 



dt 



XHel I /^Heirie + ^ CT^elCr A^^ J + XHelllttH, 

/ M 

X^ell I /^Hell^e + ^Hell^e + ^ CFmellCrNi 



= c 



XHell D. 



(A31) 



The RHS terms are, in order of appearance, Hei colli- 
sional ionizations. Hem recombinations, Hei photoioniza- 
tions (with an optional homogeneous background in paren- 
theses). Hell collisional ionizations. Hell recombinations and 
Hell photo- ionizations. In the third line we have grouped the 
terms into a creation term C and destruction terms D. 

The discrete update is done with the same formulation 



as (A30), i.e 



t+At _ X^n + CAt 
^Hell - 1+^At ' 



(A32) 



using updated values of Ar/+^*, T^+'^', x'^^^\ and x^+f * (^ 
ji and rie), and the un-updated value only of XHem- 

The only variable left is the Hem fraction. The differ- 
ential equation is 



dxY 



dt 



XHeW I /^Helirie + ^ CT^ellCr A^2 



- a^Helll <^HeIIl''^e 
-- C - XHelll D. 



(A33) 



In the third line we have as usual grouped the terms into 
creation and destruction. 

Again the update follows the same formulation. 



'^Helll 



xLiii + CAt 



l + DAt 
which is implicit in all variables. 



(A34) 
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Conservation of helium density is then enforced, i.e. that 

XHel + XHell + XHelll = 1, (A35) 

by lowering the largest of these fractions accordingly (in the 
case of XHei being the largest there is no update). 

The 10% rule is not applied to the helium fractions. 
Instead, the final 10% check is done on the electron density, 
which is retrieved from all the ionization fractions with 



Tie — XHliriH + (xHell + 2XHeIIl)nHe. 



(A36) 



If 



t+/^t 



>0.1, 



(A37) 



cool_step is restarted with half the time-step length. 



|(v) I Time-step check 

All the variables have been updated, from Ut to U^^^ , and 
the 10% rule is not violated over the time- step. The time- 
step length just used may have been unneccessarily short, 
and if so, there is a large probability that it is also unnecces- 
sarily short for the next call to cool_step (to fill the total 

AtRT). 

Therefore a final time-step check is made before finish- 
ing up, of how close we were to breaking the 10% rule. If the 
maximally changed variable in Ut has changed by less than 
5%, i.e. if 



U'^^^' 



-U\ 



U^ 



< 0.05, 



(A38) 



then the next time- step length At is twice the one just used. 
Note that this is on a cell- by-cell basis, and the next time- 
step for each cell is only stored during the subcycling and 
lost at the end of each AIrt cycle. At the beginning of each 
cell-cycle the first guess at a time-step is always AIrt- If 
this is too large for the 10% rule to be obeyed, successive 
calls to cool_step will quickly fix that by halving the time- 
step until the rule is no longer broken, and only then will 
cool.step start to return updated values oiUr- 



APPENDIX B: STELLAR UV EMISSION AND 
DERIVED PHOTON ATTRIBUTES 



In the photon injection step in RAMSES-RT (^3.1), the task 



is to inject photons into each grid cell corresponding to the 
luminosities of stellar particles that reside in it. Here we 
describe how we derive these luminosities from stellar energy 
distribution (SED) models, along with the photoionization 
cross sections and energies for each photon group. 

Bl Stellar luminosities 

Stellar particles in RAMSES represent stellar populations, so 
it makes sense to use SED models to infer their luminosi- 
ties. RAMSES-RT can read SED tables at startup and derive 
from them stellar luminosities for photon injection, as well 
as photon group attributes that can be updated to reflect 
the average emission from the stellar particles populating 
the simulation. 

We have hitherto used the SED model of IBruzual fcl 



Chariot ( 2003 ) (BC03), but it can be replaced with any other 
model, e.g. Starburst99 ( [Leitherer et al.|1999 ), as long as the 
file format is adjusted to match. The model should come in 
the form of spectra, J\ (r, Z) , giving emitted energy in Solar 
units per Solar mass per wavelength per second, binned by 
stellar population age and metallicity. In Fig. |Bl| we show J\ 
from BC03 and Starburst99 at solar metallicity for various 
population ages. 

Age and metallicity dependent population luminosity L, 
given in number of photons emitted per second into photon 
group i, is calculated from the SED model by 



where J^ — c/v'^ Jx(^^). 
ity is then 






r, Z)/hh' diy, 



(Bl) 



The cumulative population luminos- 



^i{r,Z) 



f 

Jo 



Li{t,Z) dt. 



(B2) 



Since both the photon injection and the calculation of pho- 
ton group attributes are done on the fly, Li (r, Z) and 
Hi (r, Z) must be evaluated as quickly as possible for given 
stellar particle ages and metallicities. Values of Li and ^^ 



are therefore only calculated from the SED spectra via (Bl ) 



and (B2) at simulation startup, and tabulated with equally- 
spaced logarithmic bins of age and metallicity, so that they 
can be retrieved with minimum computational effort via lin- 
ear interpolation, e.g. when injecting photons into cells via 
Eq.llTl 



B2 Photon group attributes 

There are three sets of global attributes for each photon 
group. These are average photon energies e^, average pho- 
toionization cross sections afj and energ y w eighted cross sec- 
tions afj^ that are defined in q2l(Eqs. 9|ll ). For an age and 
metallicity dependent reference spectrum Ji/(r, Z), these are 



e^{r,Z)^ 



of J (r, Z) 



r^^^ GT^jJi./hh' diy 



r^^^ (JviJv dv 

_E /_ r7\ _ Jl^iO ■J 



CTijir^Z) 



n^ J. du 



(B3) 
(B4) 
(B5) 



Since there are three ionizeable species in the current 
implementation of RAMSES-RT, each photon group has three 
values of a^ and three of a^. These attributes can be set as 
run parameters to reflect some typical stellar spectra, e.g. a 
blackbody or a SED. It can also be left to RAMSES-RT to set 
them on the fly to reflect the in-simulation stellar popula- 



tions, using the expressions (B3)-(B5), with the loaded SED 



spectra representing J^^ and the expressions from |Verner| 
|et al.| (|1996) for a^j (see Appendix E4 ) . Due to the averaged 



nature of the photon groups, we must however suffice to set 
the group attributes to reflect the average stellar emission 
in the simulation [^ If this option is used, the photon group 

^ Note that the on-the-fly update of photon attributes accord- 
ing to iB7l and (IbgI infers that existing photons attributes are 
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Chariot] ('2003') and (b) Starburst99 (Leit herer et al.|1999| > for solar metallicity at different 
The spectral luminosity is given in solar luminosities (3.8 x 10*^^ ergs~^) per solar mass (2 x 10^^ g) per 



Figure Bl. SED plots from (a) [Bruzual . 

stellar population ages. 

wavelength. Vertical lines mark the ionization wavelengths for Hi, Hel and Hell, which correspond to the wavelengths marking the three 

photon groups we typically use in our simulations. The Starburst99 spectra are generated with the instantaneous formation of 10^ solar 

masses and a Salpeter initial mass function. 



attributes are updated every n coarse time-steps (where n is 
an adjustable parameter) by polling all the stellar particles 
in the simulation and setting for each group i and species j, 



all stars 

all stars 



Y o-fj (r^ , Z^) m^ Li (r^ , Z^ 

N _ • 

•^ all stars 



(B6) 



(B7) 



Y ofj{r^,Z^ 



Li (r^ , Z^ 



^ all stars 

Y ^* Li (r^ , Z^ 



(B8) 



The values of each stellar particle's L^(T*,Z^), eij(T*,Z^), 
cr^(T^, Z-),) and afjir-),^ Z^) are interpolated from tables that 



are generated at startup via (Bl) and (B3)-(B5) 



Although one is free to use many photon groups to re- 
solve frequencies, it is practical to only use a handful, due 
to limitations in memory and computation. We typically 
use three photon groups in our simulations, representing Hi, 
Hei and Hell ionizing photons, as indicated by vertical lines 
in the plots of Fig. |B1| The stellar luminosities, instanta- 
neous and accumulated, average cross sections and energies 



for those groups are plotted in Fig. B2 for BC03 (top) and 
Starburst99 (bottom), as calculated via (Bl )-( |B5[ ). From the 
luminosity plots (top rows), it can be seen that the stellar 
populations emit predominantly for the first '-^ 3 — 6 Myrs 



changed, i.e. the attributes of photons that have already been 
emitted change in mid-air. 



and the luminosity drastically goes down as the most mas- 
sive stars in the population begin to expire. 



APPENDIX C: NON-EQUILIBRIUM 
THERMOCHEMISTRY TESTS 

To validate the non-equilibrium thermochemistry in 
RAMSES-RT we ran one-cell thermochemistry tests, that start 
at some initial state (temperature, ionization state, photon 
flux) and evolve over roughly a Hubble time. We are inter- 
ested here in verifying that our implementation is correct 
and error free and also in comparing equilibrium vs. non- 
equilibrium cooling - e.g. Cen & Fang (2006) report that 



the methods can produce significantly different results. We 
compare against the equilibrium thermochemistry of RAMSES 
which has been modified to use the exact same heating, cool- 
ing and interaction rates as RAMSES-RT. 

We test to see (i) whether the thermochemistry of 
RAMSES-RT is stable, i.e. if the stiffness of the equations re- 
sults in any sudden divergence or 'wiggles' in the evolution 
of the gas, (ii) whether RAMSES-RT evolves the ionization 
fractions towards the correct states predicted by the equi- 
librium solver of RAMSES, and (iii) whether the RAMSES and 
RAMSES-RT evolve the temperature towards the same final 
value. 

There are four tests: first we disable cooling and evolve 
only the ionization states of hydrogen and helium at dif- 
ferent constant temperatures in a zero UV radiation field, 
and see if we reach equilibrium ionization states (predicted 
by RAMSES). Then we turn on a constant UV radiation field 
and again see if we reach equilibrium states. Then we turn 
on cooling, and for two sets (zero, nonzero radiation field) 
see if the temperature evolution is comparable to RAMSES 
equilibrium cooling from the same initial conditions. 
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Figure B2. Hell, Hel, and Hi ionizing luminosities and photon group attributes derived from the BC03 (top panel) and Starburst99 
SED models (bottom panel), as functions of age (x-axis) and metallicity (colors). The plot columns represent the three photon groups. 
Top rows show stellar luminosity, in the number of photons that goes into each group per second per solar mass. Second rows show 
accumulated number of photons emitted. Tliird rows show the average photon energies per interaction. Bottom rows show average 
cross sections per interaction. 



CI Ionization convergence at constant 

temperature and zero ionizing photon flux 

In the first test, cooling is turned off and we check for a 
range of densities, temperatures and initial ionization states 
whether we get a convergence of the ionized fractions to- 



wards their equilibrium states, as predicted by RAMSES, as- 
suming zero flux of ionizing photons. 

Fig. |C1 1 shows the results. Each panel of 3 x 6 plots in 
the figure represents an evolution given the constant tem- 
perature written to the right of the panel, and shows how 
the ionized fractions, xhii, ^Heii and XHeiii, evolve from dif- 
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Figure CI. Ionization convergence test with constant T and zero ionizing photon flux. Coloured hues show non-equihbrium evolution 
of the ionization fractions, given constant T (right) and tih (top). Black dashed lines show the corresponding equilibrium ionization 
fractions as calculated in RAMSES. 
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ferent (color-coded) starting states Xi — xhii = a^Heiii (the 
Hell fraction always starts at zero). A black dashed line in 
each plot shows the equilibrium ionization fraction for the 
given temperature and species (which is gas density inde- 
pendent in the case of zero ionizing flux). Each column of 
plots represents a (non-evolving) hydrogen number density. 
The non-equilibrium ionization fractions always evolve 
towards the equilibrium ones, at a rate which depends on 
gas density, as expected. It can even take longer than the 
age of the Universe to reach equilibrium for the most dif- 
fuse gas (riH ^ 10~^ cm""^), which indeed is a significant 
difference from the equilibrium assumption. If we zoom in 
around the equilibrium states we find a difference between 
the calculated equilibrium state and the evolved one which 
is typically around one in ten-thousand - this simply cor- 
responds to the allowed error in the iterative equilibrium 
calculation, and can be decreased at will by reducing this 
error margin. 



C2 Ionization convergence at constant 

temperature and nonzero ionizing photon flux 

This is the same as the previous test, except now we apply a 
constant fiux of 10^ ionizing photons s~^ cm~^ through the 
cell, assuming the spectrum of a blackbody at 10^ K. 

Fig. |C2| shows the results. The black dashed lines in 
each plot show the equilibrium state which now is density 
dependent - the denser the gas the harder it is for the radi- 
ation field to battle against recombinations. Again the non- 
equilibrium ionized state always evolves towards the equil- 
brium one at a gas density dependent rate, though note that 
here it takes a maximum of ^ 10 Myr, which is much shorter 
than it can take in the zero photon fiux case. 



C3 Temperature convergence Avith zero ionizing 
photon flux 

Now cooling is turned on, and we compare the RAMSES-RT 
non-equilibrium temperature evolution with that of equilib- 
rium RAMSES (though keep in mind it has been adjusted to 
contain the exact same cooling rates as RAMSES-RT). Each of 
the 5 rows in fig. |C3| shows cooling for a range of decreasing 
initial temperatures, from top to bottom. The color-codings 
(initial ionization states) and columns (hydrogen number 
densities) are the same as before. The solid colored lines 
show non-equilibrium cooling in RAMSES-RT and the black 
dashed lines represent equilibrium cooling in RAMSES start- 
ing from the same temperature. 

Clearly the temperature evolution is quite similar be- 
tween equilibrium/non-equilibrium cooling, especially if the 
initial ionization fraction is 'correct', i.e. if it matches the 
equilibrium one at the initial temperature. 

The final temperature reached in the non-equilibrium 
case is usually a bit lower than in the equilibrium case. This 
is independent of gas density and initial temperature (as 
long as the initial temperature allows for cooling to hap- 
pen). The reason for this is that the non-equilibrium ion- 
ization evolution lags behind the instantaneous equilibrium 
one, so there is always a somewhat larger reservoir of elec- 
trons in the non-equilibrium case. Electrons are the primary 
cooling agents, and complete electron depletion completely 



stops cooling, so it makes sense that if the electrons deplete 
more slowly, cooling is more effective and can bring the gas 
to a lower final temperature. 



C4 Temperature convergence Avith nonzero 
ionizing photon flux 

This is the same as the previous test, except now we ap- 
ply a constant fiux of 10^ ionizing photons s~^cm~^, as- 
suming the spectrum of a blackbody at 10^ K. The results 
are shown in Fig. |C4[ Things are much the same as before, 
except that the non-equilibrium temperature seems to con- 
verge to a value which is much closer to the the equilibrium 
one - because of the ionizing fiux there is always a reservoir of 
electrons both in the equilibrium and non-equilibrium evo- 
lution, which makes for a much closer match in the final 
temperature. 

Although the final temperature reached is identical be- 
tween the two methods, the evolution towards that final 
temperature can be quite different, depending on the ini- 
tial ionization states. 

A zoom-in on one of the plots is shown in Fig. |C5| and 
reveals that there is very little difference between the final 
temperatures reached. The little difference there is results 
from interpolation from cooling-rate tables in RAMSES equi- 
librium cooling and it can be decreased further by increasing 
the size of these tables. 



C5 Thermochemistry tests conclusions 

The main conclusions of the one-cell thermochemistry tests 
are: 

• We always eventually reach the equilibrium ionization 
state with the non-equilibrium method... 

• ...but this can take a very long time to happen for dif- 
fuse gas, even more than a Hubble time. 

• Non-equilibrium temperature evolution of the gas is 
quite dependent on the initial ionization fraction of the gas 
at intermediate temperatures and low densities... 

• ...but in the end we reach the same or at least a very 
similar temperature as in the equilibrium case. 

• The convergence of the non-equilibrium solver towards 
the results of the equilibrium solver of RAMSES, given the 
same cooling rate expressions, suggests that our thermo- 
chemistry solver is robust and correct. 



APPENDIX D: ON MULTI-STEPPING IN THE 
AMR LEVEL HIERARCHY 

Solving hydrodynamics over an AMR grid with a multi- 
stepping approach always leaves ill-defined states at inter- 
level boundaries between the start and finish of the coarse 
level timestep. 

Hydrodynamic advection across the boundaries of a cell 
is performed in an operator-split fashion, such that the ad- 
vection is solved separately across a discretized timestep for 
each boundary. In order for the solver to be consistent, i.e. 
for the result at the end of the timestep to be independent 
of the order in which the boundaries are accounted for, the 
solver must work from the same initial cell state U for all the 
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Figure C2. Ionization convergence test with constant temperature and an ionizing photon flux of 10 s cm 
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Figure C3. Temperature convergence with zero ionizing flux. Color coded lines show different initial states of xhh and ^Heiii^ as 
indicated by the color legend at bottom right. Black dashed curves show the equilibrium evolution from RAMSES. 



10^ 
10' 

loH 



10^ 
10' 
10^ 



10^^ 

10' 

10^ 



10' 
10=^ 



10^ 

10't 

10= 



-H -H -H -H- 



-H -H -H -H- 




H -H -H -H— 



;^Si. 



H -H -H -H— 



:^ 



:r 



T 



1 10' 10' 1 10' 10' 1 10' 10' 1 10' 10' 1 10' 10' 1 10' 10' 

Time [Myr] 

Figure C4. Temperature convergence with an ionizing photon flux of 10^ s~^ cm~^. 
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Figure C5. Close-up of temperature convergence, for the UV 
inclusive test with initial temperature T ^ 10^ K and nn = 10~^ 
cm "^ 



inter-cell updates. Thus, a copy is first made of the original 
cell states involved, i.e. 



U^U, 



(Dl) 



where we can term lA the source state and lA the destination 
state. Using U as source terms for the intercell fluxes, the ad- 
vection can be solved with some computational method (e.g. 
Godunov solver for the hydrodynamics in RAMSES and an 
HLL/GLF flux function for the RT advection in RAMSES-RT), 
which performs the update on lA. To t ake a concrete exam- 
ple, each RT advection update, Eq. 22 uses lA for the update 
(the LHS term and the first RHS term), but the intercell 
fluxes are derived from lA^ i.e. T{U). Once all the updates (6 
per cell) have been collected, the cell update is made final 
by: 



U^U. 



(D2) 



In the cLmr_step(^) hierarchy in RAMSES, such copies are 
made of all I cells before the AMR recursion, and the up- 
date is made final after the recursion has returned and the 
hydro_solver has been called at the current level, i.e. advec- 
tion has been performed over the timestep over the current 
level and all finer levels. 

This allows cell states to be updated not only at the 
current level, but also (twice) in all neighbouring cells at the 
next coarser level. The coarser level update is only partial 
though, because it only refiects the intercell fiuxes across 
inter-level boundaries, and fiuxes across other boundaries 
(same level or next coarser level) will only be accounted for 
when the coarser level time-step is advanced. Until then, 
these coarser level neighbour cells have t wo g as states, lA 



Dl 



and lA. This is shown schematically in Fig. 

If RT subcycling is to be done at each AMR fine-level 
step, over the whole grid, the question is, which cell state do 
we use for the thermochemistry, i.e. the interaction between 
photons and gas, in those inter-level boundary cells? 

Choosing one but not the other leads to an obvious 
and severe inconsistency between the source and destination 
states. If thermochemistry does the update on lA^ then a gas 
element which is transported from one cell to it's neighbour 
during the following hydro transport is not thermochemi- 
cally evolved over the time-step, because it originates from 
lA. If instead the update is done on ^, a gas element which 
stays still in a cell over the following hydro transport step 



Figure Dl. Level i gas state updates via intercell fluxes also 
perform partial gas updates in neighbouring cells at level I — 1. 
The example shown corresponds to the hierarchy from Fig. [s] 
Steps 3 and 4 at the flnest level also include partial updates of 
neighbouring ^max — 1 cells, but these neighbour cell states are not 
fully updated until all the intercell fluxes are taken into account, 
which is in step 5 from Fig. Is] 



is not thermochemically evolved over the time-step. One 
might then just update both states via thermochemistry, 
i.e. apply it on each cell twice. This does not really make 
sense for these inter-level intercell boundary cells that have 
lA ^lA^ a,s lA doesn't represent a true state but is rather an 
intermediate and temporary quantity that exists between 
well-defined times. Also, it would be really non-trivial to 
implement: applying thermochemistry on each of the states 
also implies transporting the photons through two different 
states in each cell, which creates alternative time-lines for 
the radiative transfer! 



APPENDIX E: INTERACTION RATE 
COEFFICIENTS ADOPTED IN RAMSESRT 

Here we collect the rate coefficients used in RAMSES-RT for 
hydrogen and helium interactions, which are fitted functions 
taken from various sources. These are, in order of appear- 
ance, collisional ionization rates, recombination rates, cool- 
ing rates (collisional ionization, recombination, collisional 
excitation, bremsstrahlung, Compton and dielectric recom- 
bination), and photoionization cross sections. 



El Collisional ionization rate coefficients 



Those are in units of [cm*^ s ^] and are taken from 



Cen 



1992[), with temperature everywhere assumed in Kelvin: 



/3hi(T) = 5.85 X 10"^^ Vt {l + 



Puei{T) = 2.38 X 10"'' Vt 1 + 



/3Hen(T) = 5.68 X 10-'' VT 1 + 




-157 809. 1/T 



-285 335. 4/T 



-631 515/T 
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E2 Recombination rate coefficients 



These are all taken from Hui & Gnedin ( 1997). For readabil- 



ity, we use the following unitless functions: 

315 614 K 
T 

570 670 K 
T 

1 263 030 K 



Ahi(T) 

AHel(T) 
AHen(T) = 



The coefficients are as follows, all in units of [cm*^ s ^]: 
a^,,{T) = 1.269 X 10" 



[1 + (Ah,/0.522)0-47]i.' 



A rrTi\ Q s/ in-14 a 0.654 
ttHellU) = 3 X 10 AHei 

«Hem(T) = 2.538 xlO-'" 



Al.E 



[l + (AHen/0.522)0-47]1.923 



auu{T) = 2.753 X 10" 



agen(T) = 1.26 X 10-^" A°"^ 



«gem(T) = 5.506 X 10- 



[1 + (Ahi/2.74)0-407]2 

0.75 
Hel 



A 



1.5 
Hen 



[l + (AHen/2.74)0-407]2.242 



CHei(T) = 9.38 X 10-'' Vt 1 + 



E3 Cooling rate coefficients 

The temperature used in these coefficients is assumed ev- 
erywhere in Kelvin. Collisional ionization cooling rate coef- 
ficients [ergcm^s"^] ( |Cen||l992| ): 

Cm(T) = 1.27 X 10-'^ Vt ( 1 + ^ /^^ ^-''' '"'-'^^ 



-285 335. 4/T 



-631 515/T 



Case A and B recombination cooling rate coefficients 
[ergcm^s-^] ( |Hui fc Gnedin|l997| ): 

A 1.965 

,Ut) = 1.778 X 10- T [,^(,^,/o;41)o.50.pa.r 

men{T) = ks TOHell = fcfi T 3 X 10"" Anel^** 



(^Hen(T) = 4.95 X 10"^^ Vt 1 + 



VHemiT) = 8 X 1.778 X 10-"' T 



r]g,{T) = 3.435 x lO"''" T 



Ai 



[1 + (AHe„/0.541)0-502]2.697 



[1 + (Ah,/2.25)0-376]3.: 



•nieuiT) = ks Taieii = ka T 1.26 x 10 " AhIi^ 



VHemiT) = 8 X 3.435 x 10"^° T 



A 



[1 + (AHe„/2.25)0-376]3.72 

Collisional excitation cooling rate coefficients [ergcm'^s"^] 



( Cen|1992| ): 



^Hi(T) = 7.5x 10-'^ 1 + 



-118 348/T 



^Hen(T) = 5.54 X 10-'^ 7^-0.397 j ^ ^ 



-473 638/T 



Table E2. Phot oionizat ion energies and corresponding frequen- 
cies 

Ion species Cion ^ion 

Hi chi = 13.60 eV u^n = 3.288 lO^^ s-^ 

Hei CHei = 24.59 eV i/Hei = 5.946 lO^^ s-^ 



Hen 



CHe 



: 54.42 eV uue 



1.316 10 



16 



Bremsstrahlung cooling rate coefficients [ergcm^s ^] (jOs- 
terbrock & Ferland||2006|: 



Oiiu{T) = 1.42 X 10"^^VT 
(9Hen(T) = 1.42 X 10-'^Vt 
6'Heni(T) =4x 1.42 X IQ-^^Vt 



Compton cooling/heating rate coefficient [ergs~^] (Haiman, 
Thoul fc Loeb|1996 ) , with a the cosmological expansion fac- 
tor and Ty = 2.121 j a K the temperature of the cosmic back- 
ground radiation.: 



w(T,a) = 1.017 X 10"^^ 

Dielectronic recombination 
[ergcm^s"^] ( |Black|l98l| ): 



2.727 



2.727 



cooling rate coefficient 



CJHell(T) = 1.24 X 10" 



E4 Cross sections 



-470 000/T ^ ^ Q 3g- 



94 000/T 



Expressions for frequency dependent photoionization Hi—, 
Hell— and Helll— cross sections are used in RAMSES-RT to 
derive photon group attributes from stellar populations (Ap- 
pendix pi) . These expressions a re taken from |Verner et al.| 



(1996) (via Hui h Gnedin|l997 ) and are given in [cm ] as a 



function of photon energy e by 



cr(e) = cro \(x- \f ^ yV\ 



y 



(1 + y/yjy~ay 



;, if e ^ Cion, (El) 



(and cm otherwise), where 

e 

X = — 
eo 

and 



-?/o, 



y= ^Jx'^ + vl. 



and the fitting parameters ctq, eo, Vw^ P, Va, yo, and yi are 
given in Table |E1| The ionization energies Cion and corre- 
sponding frequencies Uion are given in Table |E2] 
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Table El. Photoionization cross section parameters - see Eg. |E1| 



Ion species 


eo [eV] 


CTQ [cm^] 


P 


Va 


Vw 


2/0 


2/1 


Hi 


0.4298 


5.475 10-14 


2.963 


32.88 











Hei 


0.1361 


9.492 10-16 


3.188 


1.469 


2.039 


0.4434 


2.136 


Hen 


1.720 


1.369 10-14 


2.963 


32.88 
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